Extended abstract

L'un des challenges majeurs de la cosmologie moderne réside en la nature même de la matière et de l'énergie noire. La matière noire peut être directement tracée à travers son effet gravitationnel sur les formes des galaxies. La mission Euclid de l'Agence Spatiale Européenne fournira précisément des données à cette fin. L'exploitation de telles données requiert une modélisation précise de la Fonction d'Étalement du Point (FEP) de l'instrument d'observation. En première approximation, les images par l'instrument des étoiles isolées donnent une mesure de sa FEP. Toutefois, ces images sont toujours bruitées, potentiellement sous-échantillonnées et la FEP sous-jacente varie dans le champ de l'instrument, mais également dans le temps et en longueur d'onde. L'objectif de cette thèse est précisément de proposer des méthodes d'estimation de FEP, à partir d'images d'étoiles, prenant en compte les différents facteurs mentionnés ci-avant, ce en s'appuyant sur des avancées méthodologiques récentes en traitement du signal.

Dans un premier temps nous avons développé une méthode permettant d'estimer à une résolution suffisante, une FEP donnée, à partir de plusieurs mesures bruitées sous-échantillonnées de la même FEP. Nous avons montré qu'en utilisant des contraintes de parcimonie et de positivité, l'on obtenait une bien meilleure restauration de la forme de la FEP, comparativement avec des méthodes largement utilisées en astronomie et en particulier à faibles rapports signal à bruit.

Toutefois, du fait de la variabilité spatiale et temporelle de la FEP, l'on ne dispose en général pas de plusieurs mesures de la même FEP. Ainsi, dans une deuxième étape, nous avons considéré le problème d'estimer conjointement un ensemble de FEPs différentes à partir d'une seule mesure potentiellement sous-échantillonnée de chacune de ces FEPs. Cela a abouti à une méthode de réduction de dimension et de super-résolution, qui en plus des contraintes de parcimonie et de positivité, exploite la compressibilité globale du champ de FEPs. Une fois de plus, nous avons obtenu un gain significatif en terme de précision sur les formes de PSFs, par rapport à des méthodes existantes.

Comme mentionné plus haut la FEP varie à travers le champ du télescope. L'objectif final étant de réaliser des mesures des formes de galaxies corrigées de l'effet de la FEP, il est nécessaire d'interpoler spatialement le champ de FEPs. En d'autres termes, l'on souhaite pouvoir estimer la FEP à une position arbitraire dans le champ de l'instrument connaissant les FEPs aux positions des étoiles isolées. Cela revient à estimer une fonction qui à tout vecteur (x, y) du plan focal de l'instrument associe un vecteur sur la variété sous-jacente aux FEPs. Pour cela, nous proposons une méthode d'interpolation qui réalise dans un premier temps une réduction de dimension locale des FEPs connues, dans la métrique de Wasserstein 2. A supposer que cela "déplie" localement la variété sous-jacente, les coordonnées des FEPs inconnues dans la représentation en dimension réduite peuvent être estimées en utilisant des méthodes d'interpolation classiques. Cela permet ensuite de calculer les coordonnées barycentriques des FEPs inconnues dans la métrique de Wasserstein 2. Les FEPs à interpoler sont finalement calculées comme barycentres des FEPs connues dans la métrique de Wasserstein 2. Nous avons obtenu par cette approche des résultats notables en termes de précision sur la forme et les valeurs de pixels des PSFs interpolées. Elle est toutefois coûteuse en termes de calculs.

Cette thèse ouvre plusieurs perspectives dont nous donnons quelques exemples ci-dessous :

-extensions des méthodes proposées pour prendre en compte la dépendance en longueur d'onde de la FEP ;

-couplage de la restauration du champ de FEP au problème de déconvolution des galaxies dans le même champ ;

-intégration davantage d'a priori physiques sur la FEP.

Notations and conventions

Conventions

All the vectors are treated as columns vectors, unless it is explicitly mentioned otherwise.

• x : a scalar.

• x : a column vector or a scalar when explicitly indicated or a manifold's point when explicitly indicated.

• x[i] : i th coefficient of x.

• X : a matrix or a random variable.

• X[i, j] : the value of the entry (i, j) of X.

• X[i, :] : the i th line of X treated as a line vector.

• X[:, j] : the j th column of X treated as a column vector.

Operators

• supp(x) : support of the vector x (set of non-zero coefficients of x).

• card(E) : cardinality of the set E (the number of elements in E).

• x p : with p > 0, p norm of x, defined as x p = p i |x[i]| p . The pseudo-norm 0 is defined as x 0 = card(supp(x)) and the infinity norm is defined as x ∞ = max({|x[i]|, i}).

• X p : with p > 0, matrix norm defined as X p = p i,j |X i,j | p (Frobenius norm for p = 2). The matrix 0 pseudo-norm and infinity norm are similarly defined.

• |||X||| : subordinate matrix norm to a vector norm • , defined as |||X||| = max v =0 Xv v .

• : element-wise matrix multiplication (Hadamard's product).

• X T : transpose of X.

• X H : hermitian transpose of X.

• •, • : scalar product.

• x : the highest integer smaller than x.

x Conventions

• rank(M) : dim(span((M[:, i]) 1≤i≤p )), where p is the number of columns of the matrix M.

• Trace(M) = n i=1 M[i, i], for M ∈ M n (R).

Notations

• N : the set of non-negative integer numbers.

• R : the set of real numbers.

• R + : the set of non-negative real numbers.

• C : the set of complex numbers.

• M np (K) : the set of n × p real or complex matrices if K = R or C respectively.

• M n (K) : the set of n×n real or complex matrices if K = R or C respectively.

• I n ∈ R n×n : identity matrix.

• L 1 (R) : the set of integrable functions on R.

• L 2 (R) : the set of square integrable functions on R.

• dim(E) : dimension of the vector space E.

• E ⊥ : the orthogonal of a vector space E.

• Span((u i ) 1≤i≤n ) : the set { n i=1 a i u i , (a i ) 1≤i≤n ∈ R n }. • Y ∼ N (m, V) : Y is a random vector that follows a gaussian multivariate distribution of mean m and covariance matrix V.

• n, m , (n, m) ∈ N 2 , n < m : the set (n + i) 0≤i≤m-n

List of abbreviations

Introduction

Sommaire

This thesis aims at proposing new methods and algorithms for astronomical images restoration, taking advantage of recent advances in applied mathematics and more specifically in calculus of variations ; a focus is made on Point Spread Function (PSF) estimation. In the following, we precise the scientific framework of our studies, underlining the challenges from a data processing perspective. Then, we give an overview of our mathematical environment. We end this introductory chapter with a description of the whole manuscript structure.

A Context

A.1 Weak gravitational lensing

The General Relativity theory can be considered as the cornerstone of Cosmology [Einstein 1915]. Over a century, it has allowed, conjointly with spectacular progress in instrumentation, an increasingly accurate understanding of the Universe general dynamic and geometry. Various cosmological models have been derived in this framework. The most widely studied and accepted is the Big Bang Cosmological model which lies on strong observational pillars such as the Cosmic Microwave Background (CMB) [Penzias & Wilson 1965] and the Universe expansion [Riess et al. 1998].

However, in order to be in perfect agreement with observations, this model has to postulate the existence of unknown and hardly observable entities, and most importantly, make them the dominant part of the Universe in terms of mass and energy. Specifically, the Universe would roughly consist of 25% percent of "Dark Matter" and 70% of "Dark Energy", the ordinary matter being confined into a few percents (see Fig. The notion of "Dark Matter" was introduced by the astrophysicist Fritz Zwicky in 1937 [Zwicky 1937]. Measuring galaxies speeds in the Coma cluster, he noticed that peripheral galaxies were excessively fast : given the observable mass of the cluster, those galaxies should have escaped its gravitational potential. He suggested as an explanation that the cluster must have an important "hidden" mass -consisting of dark matter -so that it can coherently maintain fast rotating galaxies within its gravitational potential. The Dark Energy on the other hand is directly related to the Universe global dynamic. Whereas the idea of an expanding Universe is admitted since Edwin Hubble observations in 1929, since 1998, several measurements have enlightened the acceleration of this expansion. The Dark energy presumably drives this acceleration by acting as a negative pressure.

Because of their larger proportions with respect to ordinary matter, these dark components shape the Universe at large spatio-temporal scales. Therefore, getting a better knowledge of their nature and distribution has become a major challenge for modern cosmology.

The descriptive "dark" comes from the fact that these entities can not be directly observed. However, in theory it is possible to observe and quantify the gravitational effect of the dark matter. Indeed, although it weakly or does not interact with baryonic matter, it is massive and it should bend the space-time according to Einstein equations and affect the way light travels.

Thus, a hypothetical perfectly circular galaxy would appear slightly elliptical and brighter to a hypothetical distant observer because of the integrated deflections due to dark matter, along the line of sight. This effect is known as the Weak Gravitational Lensing and is illustrated in In order to quantify this effect, one needs to be able to measure a galaxy shape from an image, which gives in practice the "lensed" galaxies shape. Since the "unlensed" galaxies are inaccessible, one needs some prior knowledge on their shapes distribution. According to the Cosmological Principle, the Universe is isotropic and homogeneous at large scales. From this principle, one infers that given a set of vectors γ i in R 2 describing the elongations and orientations of some unlensed galaxies X i , for a sufficiently large number of galaxies, one should have 1 N N i=1 γ i ≈ 0. However, the weak lensing (WL) causes an intrinsic bias in the galaxy shapes distribution which is directly related to the dark matter mass distribution and which can be estimated from the observed galaxies. The dark matter distribution in turn, can be related to the dark energy density in the Big Bang cosmological model [START_REF] Dodelson | Modern cosmology[END_REF]]. The Euclid mission, a spatial survey of the European Spatial Agency (ESA), has been designed to provide data for such a purpose, as we will see in the next section.

A.2 ESA Euclid project

Euclid is a Medium-class ESA survey mission [Laureijs et al. 2011]. It is part of the ESA's program "Cosmic Vision". Euclid primary goal is to learn about dark matter and dark energy nature. The Euclid telescope, which is space-based, will be launched on a Soyuz rocket from Europe's spaceport in Kourou, in Guiana ; the launch is currently planed for December 2020. Quantifying WL requires measuring accurately shapes over a large set of galaxies. To that end, the Euclid telescope -provides a large field-of-view (fov) ; the embedded visible and nearinfrared instrument have a fov of roughly 0.54deg 2 which is two orders of magnitude higher than Hubble Space Telescope (HST) camera's ;

-is extremely sensitive with a lower magnitude bound of AB mag 24.5 (this corresponds to a spectral luminosity 17 order of magnitude lower than the sun's) ; this enables the observation of distant galaxies which is essential, since the further is a galaxy, the more it is potentially lensed ;

-is planned to orbit around the Sun-Earth Lagrange point 2 (see Fig.

I.
3) ; consequently it will be free of atmospheric turbulence (which is true for any space-based telescope), as well as Earth and Moon magnetic fields disturbances ; moreover, the payload will neither be occulted nor illuminated by the Sun, the Earth or the Moon during the survey ; finally, there is no significant variation of gravity along the selected orbit ; all these features will yield an exceptionally homogeneous final data-set ;

-is extremely stable in terms of structural vibrations, with a specified pointing jitter in average smaller than the pixels size. The survey active phase will last 6 years, providing at the end a legacy of more than 1 billion of well resolved galaxies images, positions and speeds, distributed over 15000deg 2 of sky without the milky way and roughly 10 billions light years in depth, thanks to its extreme sensitivity. Henceforth, Euclid will show the Universe at a spatio-temporal scale that has never been reached before, which makes this mission exciting and promising for Cosmology.

Regardless of how clean will Euclid data presumably be, it is crucial however to have an accurate knowledge of the instrumental effects on the galaxies shapes, in order to exploit them properly. Indeed, the image forming process itself induces shape distortions that might be mistaken for the Weak Lensing (WL) effect. This is one of the biggest challenge in WL measurement, as we detail in the next section.

A.3 Data and instrumental effects

In this section, we review the image forming process for an Euclid-like telescope, precising the distortions introduced with respect to the "real world objects". Despite galaxies are 3D objects, our conceptual real world objects are 2D light distributions contained in a plane which is perpendicular to the line of sight as illustrated in Fig. I.4 ; we do not consider the geometrical distortion associated with this planar projection. projection with respect to the real world 3D objects is not considered.

A.3.1 Optics related effects

The optical instrument warping effect is largely related to the phenomenon of diffraction, which characterizes a propagating wave behavior whenever it encounters an obstacle. Therefore, we make a brief recall of diffraction theory applied to optics thereafter. The formulae are extracted from [Goodman 2005].

According to the Huygens-Fresnel principle, when a wave propagates, each point of the space reached by the wave becomes in its turn the source of a spherical wave. This is illustrated in Fig. I.5.

These secondary waves can interfere : according to their phases, they can add up in amplitudes or attenuate one another. In particular, when the wave encounters an obstacle, as one may expect intuitively, the secondary sources physical properties change in the vicinity of the boundaries so that the overall transmitted wave might significantly depart from the incident one ; this is especially true when the size of the aperture (see Fig. The first rigorous description of the diffraction by a planar aperture in terms of superposition of secondary spherical waves is due to the physicist Gustav Kirchhoff (1882). He proposed a mathematical framework built upon the Greens theorem and derived the diffracted wave making some simplifying assumptions known as the Kirchhoff boundary conditions. These assumptions were later proved to be inconsistent by Poincaré (1892) and eliminated by Sommerfeld (1894) who proposed the so-called Rayleigh-Sommerfeld diffraction theory. The latter characterizes the diffracted wave in a parallel observation plane located at normal distance z from the aperture (see Fig. where p 1 = (x, y, z), p 0 = (η, ξ, 0), r 01 = p 1 -p 0 2 , U (p 0 ) is the incident wave amplitude and phase at the point p 0 of the aperture, λ is its wavelength, k = 2π λ and Σ is the aperture. Two interesting approximations of the Rayleigh-Sommerfeld diffraction are later proposed. The first one is the Fresnel approximation. It relies on the following approximation :

r 01 ≈ z(1 + 1 2 ( x -η z ) 2 + 1 2 ( y -ξ z ) 2 ). (I.2)
Taking r 2 01 ≈ z 2 in I.1 and replacing the r 01 appearing in the exponent with I.2, one can express the diffracted field as U (p 1 ) = e jkz jλz e j k 2z (x 2 +y 2 ) R 2

(i Σ (p 0 )U (p 0 )e j k 2z (η 2 +ξ 2 ) )e -j 2π λz (xη+yξ) dηdξ, (I.3)

where we can identify the Fourier transform of the incident wave delimited by the aperture and multiplied by a quadratic phase term. The regime in which this approximation holds is known as the region of Fresnel diffraction, or the near field of the aperture, since it accounts for the curvature of the diffracted wavefront through the quadratic phase term. It yields accurate results when the major contributors in the integral I.1 are the terms for which x ≈ η and y ≈ ξ.

Besides, the Fraunhoffer approximation additionally assumes that z kmax Σ (η 2 + ξ 2 ) 2 , (I.4) which implies that the quadratic phase term in the integral in I.3 is approximately equal to 1 everywhere in Σ. Therefore, the diffracted field expression simplifies to

U (p 1 ) = e jkz jλz e j k 2z (x 2 +y 2 ) R 2 i Σ (p 0 )U (p 0 )e -j 2π λz (xη+yξ) dηdξ. (I.5)
The diffracted wave is given by the Fourier transform of the aperture limited incident wave up to a multiplicative term. Now we consider a general linear optical system that we treat as a black box as represented on Fig.I.7.

The optical system maps an incident scalar wave U i into an output wave given by the so-called superposition integral : where h o (x, y, u, v) is the field's value at the location (x, y) of the exit pupil that one would get from the optical system as a response to an unitary impulse located at (u, v) in the entrance pupil. h o is the impulse response of the system. The system is stationary if for a couple of locations {(u, v), (x, y)}, h o (x, y, u, v) = h o (u-x, v-y) ; in words, shifting an impulse in the object space translates into a shifting of the output wave, with no changing in its form. For stationary systems, the right hand side in Eq. I.6 becomes a convolution. Real imaging systems are never perfectly stationary and this non-stationarity motivates an important fraction of this thesis work. The concepts of entrance and exit pupils generalized the notion of aperture used in deriving the Rayleigh-Sommerfeld and Fraunhoffer diffraction approximations. Firstly introduced by Ernst Abbe, the entrance pupil is the image of the most severely limiting aperture when viewed from the object space ; conversely, the exit pupil is the image of the same limiting aperture viewed from the image space. Regardless of the optical system's complexity, the wave that reaches the image's plane can be described based on the diffraction effects that would result from the exit pupil. We assume that the light propagation between the entrance and exit pupil plans can be well described by geometrical-optics. The system is said to be diffraction-limited. The ideal geometrical-optics predicted image of the incident wave can be defined as

U o (x, y) = R 2 U i (u, v)h o (x,
U g (u , v ) = 1 |M | U i ( u M , v M ), (I.7)
where M is a magnification factor. Under these hypothesis, the image wave is given by

U o (x, y) = R 2 U g (u , v )h o (x -u , y -v )du dv , (I.8)
where the impulse response is calculated within the Fraunhoffer approximation as

h o (u, v) = A λz i R 2 P (x, y)e -j 2π
λz i (xu+yv) dxdy. (I.9)

A is a constant amplitude. z i is the distance from the exit pupil plan to the image plan. P is a binary mask defining the exit pupil. It is often referred to as the pupil function. In most imaging system, the photo-current generated by the sensors is proportional to the incident power density. Therefore, one can measure to the diffracted field's intensity which is defined as

I(x, y) = |U o (x, y)| 2 (I.10)
for a monochromatic light source or more generally

I(x, y) = lim 1 T T /2 -T /2 |U o (x, y, t)| 2 dt = |U o (x, y, t)| 2 t (I.11)
for polychromatic sources. This last definition is particularly relevant from a practical of view since the photo-detectors integration time are in general substantially longer than the incident wave temporal period. Assuming that the field U i has a sufficiently narrow spectrum, the impulse response h o can be considered to be wavelength independent so that Eq.I.8 is generalized as follows :

U o (x, y, t) = R 2 U g (u , v , t -τ )h o (x -u , y -v )du dv .
(I.12)

The delay τ is related to the propagation from (u , v ) to (x, y). One can justify that under the narrowband assumption, the intensity can be written as

I(x, y) = R 2 du 1 dv 1 R 2 du 2 dv 2 h o (x-u 1 , y-v 1 )h * o (x-u 2 , y-v 2 )J g (u 1 , v 1 , u 2 , v 2 ) (I.13) where J g (u 1 , v 1 , u 2 , v 2 ) = U g (u 1 , v 1 , t)U * g (u 2 , v 2 , t) t .
J g is the so-called mutual intensity and measures the spatial coherence of the field U g .

If the illumination is spatially coherent, the phases of the different impulses in the object's plane are perfectly correlated. This is nearly the case if the light is emitted by a laser. Reversely, the illumination is said to be spatially incoherent if the phases vary in an uncorrelated way in the object's plane. Most of the common light sources, including natural ones can be in first approximation considered as incoherent. In this case, the mutual coherence is simplified to

J g (u 1 , v 1 , u 2 , v 2 ) = I g (u 1 , v 1 )δ(u 1 , u 2 )δ(v 1 , v 2 ), (I.14)
where I g is the intensity associated with the field U g . From Eq.I.14 and Eq.I.13 one obtains the following important equation relating the incident and diffracted field intensities for a perfectly incoherent lighting :

I(x, y) = R 2 |h o (x -u 1 , y -v 1 )| 2 I g (u 1 , v 1 )du 1 dv 1 . (I.15)
In words, for a diffraction-limited imaging system, when the illumination is incoherent, the image intensity is obtained by convolving the ideal geometrical image intensity with the intensity impulse response |h o | 2 , which is mostly referred to as the Point Spread Function (PSF).

The Optical Transfert Function (OTF) is defined as the normalized Fourier transform of the PSF. For a diffraction-limited imaging system, it can be shown that the OTF is the normalized autocorrelation function of the pupil function P (see Eq.I.9). As such, the OTF has some interesting properties among which it is an even function. The pupil function being by definition spatially bounded, this also implies that a PSF is bandlimited.

The results stated so far were derived assuming that the optics do not introduce any phase-shifting as the incident wave goes though the system. However these phase-shiftings also termed as wavefront errors always occurs due to optical aberrations. In order to account for the wavefront errors, one shall generalize the notion of pupil function. We note the phase-shifting at a point (x, y) φ(x, y) = 2π λ W (x, y), where λ is the central wavelength of the spectrum of the incident wave. W (x, y) represents the optical path difference. Then, the generalized pupil function is defined as P(x, y) = P (x, y)e j 2π λ W (x,y) . (I.16)

Then the amplitude impulse response is given by the diffraction pattern of an aperture having a complex amplitude transmittance P. The PSF's definition remains unchanged.

It is worth noting that the pupil function P actually depends on the position of the light source position w.r.t. the instrument fov. This makes the PSF spatially variable. In general, computing the pupil function requires performing a ray tracing using a dedicated software as Zemax. Moreover, the different optical aberrations effects are more or less pronounced depending on the light source's position in the fov. Let consider the well-known coma aberration for instance. The coma is due to a variation of magnification across the fov. Indeed, when a bundle of rays coming from a single off-axis source hits a lens or a mirror, each ray is reflected toward a different point in the focal plane. The more a ray is close to the optical element edge, the further from the optical axis the reflected ray crosses the focal plane. This gives isolated star images a cometary shape, hence the name of this aberration. As illustrated in Fig. I.8, the coma varies across the fov, which is translated by a spatial dependency of the function W is Eq.I.16. Thus optical aberrations contribute Figure I.8 -Simulated coma aberration : the coma varies across the fov ; in particular, the further is the source from the fov center, the more dramatic is the distortion. Credit : nikon.com to the PSFs variations across the fov. The optical path difference function should be noted W x i ,y i to explicit its dependency to the impulse's position in the object's plane. However we omit these indexes for simplicity. We shall add that aberrations such as defocusing or tilt can vary in time due to thermomechanical constraints on the imaging system. This is especially relevant for space based telescopes. This makes the PSF ultimately time dependent. Moreover, as shown by Eq.I.9, the PSF also depends on the wavelength ; for a single channel imaging system, the image of a broadband source of light is the result of a continuous summation of ideal monochromatic images convolved with a PSF that might considerably change across the spectrum.

A detailed presentation of optical aberrations can be found in [Mahajan 1998].

A.3.2 Electronics related effects

So far, the images have been modeled as continuous intensity distributions. However for numerical images, the diffracted intensity distribution is recorded via a matrix of photo-detectors. We assume that the detectors cover 100% of the focal plane. Each detector integrates the light intensity over its surface during a certain amount of time. We note ∆ the detector matrix step and D the focal plane area covered by a reference detector. The intensity recorded by the detector (k, l) is in first approximation given by I kl = D I(k∆ + x, l∆ + y)dxdy, (I.17) up to a multiplicative constant. Therefore, the numerical image forming process can be described as a convolution of the intensity distribution I with the kernel defined as

h D (x, y) = 1 if x ∈ D,
0 otherwise, (I.18) followed with a sampling. However, this detector model is idealistic and we give below a nonexhaustive list of limitations :

-thermal noise : even in the absence of illumination, the detectors unavoidably generate a signal due to random thermal motions of electrons ; the resulting "background image" is generally modeled as a white gaussian process [Nyquist 1928] ;

-dark-current shot noise : this noise is due to random generation of electrons in the detectors and is also related to the temperature ; however unlike the thermal noise it is non-gaussian [Baer 2006] ;

-readout noise : this noise is inherent to CCD detectors ; it is related to the uncertainty on the photo-electrons total charge measurement in the CCD matrices ships due to the imperfections of the electronic components ; it is well modeled as white gaussian process [Basden et al. 2004] ;

-charge transfer inefficiency (CTI) : this effect is particularly present for space telescope ; the high energy radiations might gradually damage the CCD detectors ; this might cause delays in the electrons transfers to the CCD matrices ships which produces a trail and blurring in the final images [Massey et al. 2014] ;

-quantization : in other to be stored numerically, the measured intensities have to be represented using a finite numerical "alphabet", which implies a quantization of the measurements and produces what is known as the "quantization noise" [Widrow & Kollár 2008].

Ignoring some of the above effects, we can summarize a digital image forming process as follows :

I d = S(h D * h o * I g ) + N, (I.19)
where S is a sampling operator and N is the noise. In the following, we refer to the function h D * h o as the PSF. A natural question arising then is whether or not the sampling operation preserves the continuous light intensity information content. In other words, can the continuous function h D * h o * I g be calculated from the discrete array S(h D * h o * I g ) ?

The Nyquist-Shannon sampling theorem gives a sufficient condition. We recall its classical formulation : if f ∈ L 1 (R) and if its Fourier transform f's support is included in the interval [-B, B], then

f (x) = lim n→+∞ k=n k=-n f k 2B sinc 2πB x - k 2B . (I.20)
In words, if an integrable function is bandlimited with a maximum frequency B, then it can be "reconstructed" from an uniform sampling with a step 1 2B . This theorem naturally extends to higher dimensions. In practice, the light intensity distribution that reaches the detectors matrices is by construction spatially bounded so that the integrability condition always applies. The second condition is guaranteed by the fact that the PSF is bandlimited.

Therefore, it is possible to preserve the focal plane intensity information, up to the electronics defects, if the sampling step, which is determined by the detectors size, is sufficiently fine. As we detail in Chapter IV, an insufficient sampling rate might cause frequency aliasing, which is a frequency content distortion independent of the instrumental effects mentioned so far.

A.3.3 Other sources of images distortions

We consider external sources of distortions. In astronomy, forming images requires long exposure times because of the weak luminosity of most galaxies seen from the solar system. For instance, the nominal exposure time for the Euclid telescope is 540s [Laureijs et al. 2011].

Due to the intrinsically random nature of the light emission process, the number of photons that hits a detector during a given amount of time can be modeled as a Poisson process. The uncertainty related to the photons count is what is referred as shot noise.

During the exposure time, a space telescope is not perfectly steady which yields a convolutive motion blur which is time varying. For ground based telescopes, the atmospheric turbulence causes a blurring which is both space and time varying. The Euclid telescope will be put into orbit above the atmosphere so that the latter blur need not to be considered. The telescope pointing jitter is specified not to exceed a pixel size, which makes the motion blur moderated ; besides, it can be readily accounted for since pointing direction times series will be available.

In this work, we consider monochromatic images. The Chapters IV to VI are dedicated to estimating discretized but well sampled PSFs in the telescope's field from undersampled and noisy stars images and assuming that the detectors response is uniform and linear. The PSFs spatial variability is taken into account. The quantization noise is considered to be negligible. The shot noise which will be more or less considerable depending on the brightness of each particular observation is not treated. The CTI effect as well is not considered.

B Methods and tools

A physical model free framework for modeling a telescope's PSFs field is proposed. This framework aims at extracting the information available in a set of instrumental data at best, while making fairly general regularity assumptions on the PSFs. We give a brief overview of the key concepts and tools in the following.

B.1 Low complexity

In the Fraunhoffer approximation, the amplitude impulse response of a linear imaging system is a scaled Fourier transform of the pupil function, up to a multiplicative factor (see Eq.I.9). This makes the optical PSF a uniformly continuous function. Therefore a discrete well sampled PSF is at least a piece-wise smooth signal, especially if it includes the detectors response. This structure or regularity implies that a discrete well sampled PSF actually has less degrees of freedom than the number of pixels needed for a direct representation. It is in principle possible to represent the PSFs in a more concise manner. This notion of conciseness classically referred as sparsity has been extensively studied and used for solving inverse problems. It will be shown very useful, especially regarding noise robustness. Besides, we mentioned that the PSF varies spatially, due to the changes of the generalized pupil function across the fov (see Eq.I.16). More specifically, the modulus function P varies according to the whole optical instrument complexity and the optical path difference function W variations are related to the optical aberrations and therefore, the optics quality. In Euclid, the spatial variations of the optical PSF are mostly due to the aberrations. We can roughly distinguish two fundamentally different types of aberrations :

-the global aberrations related to the global structure and defects of the optics as well as their relative positioning ; they include effects such as defocusing, tilt, coma etc. ; they induce a smoothly varying phaseshifting across the fov and constitute the low frequencies of the function W ; these aberrations are well-modeled using Zernike polynomials [START_REF] Wang | [END_REF] ;

-the local aberrations due to polishing defects ; indeed, the small residual structures distributed over the optics surfaces diffuses light inducing random phase-shiftings ; thus these aberrations constitute the high frequencies of W and their effect is susceptible to change rapidly w.r.t. the impulse's position in the object's plane.

If we assume that the wavefront errors are mainly determined by the global aberrations, then the PSFs evolves smoothly across the fov and are strongly correlated. The latter point implies that a set of discrete well-sampled PSFs is scattered along a limited number of directions relative to their ambient space dimension. This low dimensionality assumption is crucial, especially in the perspective of a joint PSFs field super-resolution. The smooth spatial variation of the PSFs links the aforementioned low dimensional scattering of the PSFs to their distribution in the fov. Considering a finite set of PSFs, this link can be elegantly built through an undirected graph inversely weighted by the PSFs locations pairwise distances in the fov. Then, the graph's function associating each node to each PSF can be enforced to be smooth by seeking sparsity.

B.2 Manifold learning

Loosely speaking, a d dimensional manifold is a topological space that can be injectively mapped to R d locally. d represents the intrinsic dimension. This definition encompasses linear subspaces, also called linear manifolds, but also non-linear sets as we will see in the following. If a manifold is a collection of real world data, the intrinsic dimension can be interpreted as the number of degrees of freedom of the physical system that generated those data. For example, Fig. I.9 shows samples from the manifold consisting of images of a given subject face with a constant pose, for all the lighting angles and positioning possible. Although the ambient space dimension, which is the number of pixels, can be substantially large, the intrinsic dimension is limited to five : three spatial coordinates and two angular coordinates.

For a fixed wavelength, the manifold model naturally arises in our study since the optical PSF is completely and uniquely determined by two parameters which are the impulse coordinates in the object plane. As shown in Section A.3.1, there is a highly non-linear relationship between these coordinates and the PSFs, which makes the manifold perspective particularly interesting.

The Manifold Learning (ML) is a discipline at the confluence of geometry and data science which aims at uncovering information on a manifold from a set of data sampled from this manifold. For example one can estimate the manifold's intrinsic dimension, learn a low dimensional parametrization or characterize the manifold's curvature. Let recall that ultimately, we want to estimate the PSFs at galaxies locations in the telescope fov, where no PSFs measurements are available. This requires being able to perform displacement .10 shows a plot of an "S-curve" data set, which is regularly used for ML algorithms benchmarking. It is clear from this plot that the euclidean distance between two points far apart (in the euclidean sense) might strongly depart from the minimal length of the paths joining these points on the manifold.

By contrast, the euclidean distances between neighbor points is approximately equal to the geodesic distances between these points, if the manifold sampling is sufficiently dense. Thus, most manifold learning techniques build over the hypothesis that the geometry of the manifold is approximately euclidean locally and thereby are sensitive to the curvature/sampling density trade-off.

In the PSFs interpolation framework introduced, we propose to replace euclidean distance with an alternative distance based on Optimal Transport (OT) theory for characterizing the local PSFs manifold geometry. Indeed, Optimal Transport offers an elegant approach for comparing positive unitary mass distributions, robustly to non-linear warping. Moreover, from a practical point of view, the Optimal Transport is particularly appealing in a context of ML for interpolation since recipes and algorithms for computing geodesics

C Manuscript structure

This manuscript is organized as follows :

-in Chapter II, we make a formal presentation of low complexity data models, focusing on sparsity and low rankness ; we describe how these priors on the data structure may be used in solving inverse problems ;

-Chapter III is dedicated to geometrical concepts and tools useful in handling data lying on a non-linear manifold ; we give a general overview of ML, followed by an introduction to the OT and its geometrical features.

-in Chapter IV, we propose a sparsity-driven multiple frame PSF superresolution (SR) method ;

-a further step is made in Chapter V where we introduce a single frame PSFs field super-resolution ; the requirement of having several low resolution (LR) observations of the same PSF is dropped and compensated by the PSFs field regularity previously mentioned ;

-the PSFs field interpolation is treated in Chapter VI ; we build a methodology for estimating a monochromatic PSF at an arbitrary location in a telescope from a set of PSFs at known locations, using Manifold Learning ideas and Optimal Transport tools ;

Chapitre II In numerous problems in signal processing, the goal is to estimate an unknown signal from an altered observation. It is in general assumed that the signal of interest follows some sort of model or has a specific structure which knowledge can be used in the estimation process. The signal's structure can often be described using a notion of simplicity or compressibility which translates either the stability in space or time of the physical system that generated the signal at a given scale of observation, or the generating system limited number of degrees of freedom. This is illustrated in Fig.II.1. This intrinsic simplicity of the data, more often referred to as low complexity prior in the signal processing literature is at the core of this chapter. After a general introduction to linear inverse problems, we formally present some important low complexity data models and give a preview of how they can be practically used for signal restoration in a variational framework. where N is the electronic related noise, I g is a physical continuous light distribution characterizing the real world scene, K is the optical device's PSF and S is a function from I to R p realizing the sampling, I being any functional space that contains I g (for instance I = L 1 (R 2 )). While the continuous func-tions K and I g are practically inaccessible, one shall assume that they can be well approximated with piece-wise constant functions over a given regular grid in R 2 . Because of the finite extension of the focal plane, we can assume that the support of I g is comprised in [0, l] × [0, l] for some positive real l chosen in such a way to also encircle most of the PSF's energy. We define the function

Low complexity data models in inverse problems

A Linear inverse problems

h 0,0 (x, y) = 1 (l/n) 2 if (x, y) ∈ [0, l/n] × [0, l/n], 0 otherwise, (II.2)
for some positive integer n. We define the orthonormal family in

L 2 (R 2 ) (h i,j ) 0≤i,j≤n-1 by h i,j (x, y) = h 0,0 (x -il/n, y -jl/n). (II.3)
Then I g and K can be approximated with their projections on Span((h i,j ) 0≤i,j≤n-1 ), which is a n 2 dimensional subspace :

Îg = 0≤i,j≤n-1 a ij h i,j , (II.4) K ≈ 0≤i,j≤n-1 b ij h i,j .
(II.5) (a ij ) 0≤i,j≤n-1 and (b ij ) 0≤i,j≤n-1 are digitally tractable representations of I g and K respectively. Moreover, for (r, s) ∈ 1, 2n -1 2 ,

( Îg * K)(rl/n, sl/n) = (i 1 ,i 2 )∈[0,n-1] 2 / i 1 +i 2 =r-1 (j 1 ,j 2 )∈[0,n-1] 2 / j 1 +j 2 =s-1 a i 1 j 1 b i 2 j 2 , (II.6)
where one recognize a discrete convolution of A = (a ij ) 0≤i,j≤n-1 and B = (b ij ) 0≤i,j≤n-1 on the right hand side. n is a parameter set by the practitioner according to the resolution needed and typically, n 2 ≥ p. For a well chosen n, we can write S( K * Îg ) = D(A * B), where D is a discrete downsampling operator.

Finally, the forward observation model becomes

Im d = D(A * B) + N, (II.7)
where N accounts for the electronic related noise and the modeling error. The task of estimating A given Im d and knowing B is a linear inverse problem because of the linearity of the degradation operator D(. * B) and because one is going backward from the observation to its cause. Linear inverse problems occur in a wide range of applications :

-in geophysics where knowledge of Earth deep structure is inferred from Earth surface seismograms [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF]] ;

-in acoustics where one can characterize the propagation medium heterogeneity, namely the obstacles, from measured waves [START_REF] Colton | [END_REF] ;

-in quantum mechanics where nuclear forces can be characterized based on observable scattered radiations in particles colliders [Chadan & Sabatier 2011] ;

-in medical imaging where one is to construct an image of something which is in the human body from ideally non-invasive and nondestructive and hence indirect measurements ; there exists several medical imaging modalities among which one of the most important is the Computed Tomography [Herman & Sabatier 1987] ;

-in astrophysics, as presented in [Starck 2016] and the references therein.

A linear inverse problem can always be written into the following generic form : In all this manuscript, • ≡ +, and the noise is said to be additive. However, coherent imaging systems such as Synthetic Aperture Radars or medical ultrasounds are predominantly affected by a multiplicative granular noise known as speckle [Gagnon & Jouan 1997].

y = Mx • b, (II.8) where -y ∈ R p is the observed signal, -x ∈ R m is the unaltered signal, -M ∈ M pm (R)

A.2 Ill-posedness and ill-conditioning

As previously stated, solving the inverse problem II.8 is the task of estimating x from the observation y. Without any knowledge on x, the accuracy of this estimation is tightly related to how much information is lost along the cause-effect sequence modeled by the matrix M (see the interesting discussion in [Bertero & Boccacci 1998]).

If M is invertible, then the equation y = Mx (i) admits a solution which (ii) is uniquely determined by xb = M -1 y. Moreover, (iii) xb continuously depends on y. These conditions of existence, uniqueness and continuity makes the problem II.8 well-posed in the sense of Hadamard [Hadamard 1902, Hadamard 1923]. The condition (iii) is meant to guarantee that one can get a physically meaningful solution from imperfect observations ; indeed, a sufficiently moderated noise implies a small departure of xb from the ideal solution corresponding to b = 0. However, M is in general non-invertible, and different situations can occur. If m > rank(M) then this equation admits multiple solutions, whenever a solution exists. This is the case when m > p, for instance if M realizes a masking of the underlying signal. If p > rank(M), then the equation y = Mx doesn't necessarily admit a solution. The problem is said to be ill-posed. In general one is more concerned with the solutions multiplicity than with the strict invertibility of M ; indeed, it is not suitable to actually solve the equation y = Mx because of the noise. Besides, even if M is perfectly invertible, estimating x might be tricky from a numerical point of view, depending on M's conditioning number, which is defined thereafter. The condition number of M is defined as

κ(M) = |||M||| M -1 , (II.9)
where |||•||| is a subordinate matrix norm to some vector norm • in R m . Assuming that y 1 = Mx 1 and y 2 = Mx 2 , the following important relation holds :

x 2 -x 1 x 2 = κ(M) y 2 -y 1 y 2 .
(II.10)

If κ(M) takes a large value, then a small change in the observation space can translate into a large change in the undistorted signal space, and the matrix M is then said to be ill-conditionned. This is typically the case in deconvolution problems where M represents a convolution with a low-pass filter. As illustrated in Fig. II.2, even with a very small amount of noise, the signal estimated by direct inversion of M can dramatically depart from the truth, even if strictly speaking, the solution of y = Mx continuously depends on the observation. Ill-conditioning besides implies an amplified sensitivity to round-off errors due to the finite machine numerical precision.

Most inverse problems are ill-posed and/or ill-conditioned. Both cases are illustrated in Fig.

II.3.
This requires the practioner to regularize the problem which is making use of prior knowledge of the signal of interest in order to narrow down the domain of possible solutions and stabilize the inversion process. We present some classical regularizers in the next section.

A.3 Regularization

A.3.1 Approximated problem

As precised in the previous section, given an observed signal y, it is not suitable to exactly solve the equation y = Mx because of noise and modeling imperfection. However one wants y and Mx to be close in some convenient sense that depends in general on the noise statistics. Precisely, the departure is measured as C(y, Mx), where C is a positive cost function which verifies

∀(u, v) ∈ R p × R p C(u, v) = 0 ⇒ u = v.
(II.11) Thus a perfect solution can be recovered by minimizing J(x) = C(Mx, y) when the problem is well-posed.

Example A.1. -In numerous applications, the noise is assumed to be a realization of a gaussian random vector B ∼ N (0, σ n I p ). Thus, the observation y is a realization of a gaussian random vector Y ∼ N (Mx, σ n I p ). Under this probabilistic model, the likelihood of x immediately gives one the popular quadratic cost C(y, Mx) = 1 2 y -Mx 2 2 . More generally, if B ∼ N (0, V) for an arbitrary positive definite correlation matrix V, then the log-likelihood derived cost function is

C(y, Mx) = 1 2 (y -Mx) T V -1 (y -Mx).
-If additionally the observation contains spurious artefacts (for example cosmic rays impact print in astronomical images), it is suitable to require a weaker data fidelity on the largest observation entry. This leads one to outliers-robust cost functions such Huber loss [Chen et al. 2014] :

C(y, Mx) = 1 2 y -Mx 2 2 if y -Mx 2 2 ≤ β, β y -Mx 1 -β 2 2 otherwise.
(II.12)

-Often, digital images can be corrupted by an additive impulsive noise due to malfunctioning optical sensors, transmission errors, etc. Such noise can be model as a realization of a laplacian random vector, yielding an Mx][i] y[i] ) (see for example [Dupé et al. 2011]) ; this cost is in general suitable when handling positive data.

We use a quadratic cost in the chapters IV and V. However astronomical images are actually corrupted with a mixed Poisson-Gaussian noise. Still, the corresponding likelihood-derived cost can be readily integrated into the proposed restoration schemes [Jezierska et al. 2013]. Note that even if the function C is a strictly convex function, the cost C(y, Mx) might have multiple global minimum in x if the problem is ill-posed. Since then, the inversion task consists in selecting a solution in the vicinity of one of a minimizer of the cost, based on prior information on the signal to be estimated. This is detailed in the following section.

A.3.2 Examples of regularizers

Tikhonov regularization One of the simplest information that can be used in an inverse problem is that the restored signal's energy, measured by its squared 2 norm, can not be too large. This can be done by solving the following optimization problem :

min x 1 2 y -Mx 2 2 + λ x 2 2 , (II.13)
where the regularization parameter λ balances the strength of the constraint on the solution's energy against the accuracy of the data fit ; intuitively, one can anticipate that the practical choice of λ will be tight to the noise level. This is a particular instance of the widely used Tikhonov regularization method. Besides, if M represents a circular convolution, it can be shown that the solution of II.13 is the optimal estimate of x in the sense of Wiener [Levinson 1946] for the particular choice λ = SNR -1 . Problem II.13 can be turned into the more generic following form, (II.14) where the matrix Γ allows for making use of diverse a priori. This regularization method became particularly popular due to its application to integral equations by the mathematician and geophysicist Andrey Nikolayevich Tikhonov [Tikhonov & Arsenin 1977]. Typical choices for Γ are finite difference approximations of gradient and laplacian operators which in both cases allows for limiting high frequencies energy in the restoration process. Problem II.14 can be interpreted in a Bayesian framework. Indeed, assuming that the unobserved signal x is a realization of a random vector X verifying ΓX ≡ N (0, σ s I p ), (II.15) then problem II.14 is equivalent to a Maximum a Posteriori (MAP) estimation of x, for λ = σ 2 n σ 2 s . We recall the probabilistic model of the noise : B ∼ N (0, σ n I p ). If Γ is a discrete approximation of a gradient operator, then II.15 reads as a brownian motion-like model on the latent signal [Karl 2005], which clearly misses local correlations inherent to natural images for instance. As anticipated in the case Γ = I p , Tikhonov regularization yields a signal estimate which depends linearly upon the observation. This is convenient since the uncertainty on the estimate is straightforwardly characterized, which is particularly suitable in cosmology for large scale statistical studies [Bobin et al. 2015]. However, high frequency informative content and noise are identically attenuated (or preserved), yielding overly smoothed (or noisy) solutions. This limitation motivated the study of non quadratic regularizers. We consider the general formulation of the estimation problem

min x 1 2 y -Mx 2 2 + λ Γx 2 2 ,
min x 1 2 y -Mx 2 2 + λR(x), (II.16)
Tikhonov regularization being a particular instance.

Some non quadratic regularizers One of the most widely used non quadratic regularizer is certainly the entropy defined as

H(x) = - p i=1 x[i] log(x[i]), (II.17)
where x is assumed to be a discrete probability distribution. The choice R(x) = -H(x) is particularly sound when the only prior available is the positivity (see [START_REF] Dominikus | Restoration of degraded images with maximum entropy[END_REF]] and the references therein). Indeed, maximizing the entropy yields a solution with less structural information or equivalently it enforces pixels values similarity. According to Shannon's interpretation of the entropy from the information theory perspective, this yields the most "natural" solution, whenever the positivity prior applies. Although this constraint implies a smoothness of some sort, it doesn't limit the frequency band of the signal to be restored, making sharp details restoration in principle possible, while smoothing the noise out [Frieden 1972]. A thorough presentation of entropy based regularization for image processing can be found in [Starck & Murtagh 2006]. However, by construction, this approach doesn't specifically account for natural signal or images properties, among which the presence of discontinuities. This leads us to the important total variation regularization. Firstly used for image denoising in the foundationnal work by Rudin, Osher and Fatemi [Rudin et al. 1992], total variation based restoration has recently known a renewed interest in the signal processing community due to theoretical and algorithmic advances ; see [Chan et al. 2006] and [START_REF] Chambolle | [END_REF] for comprehensive reviews. Total variation has been defined and theoretically analyzed mostly in a continuous setting ; however the following discrete approximation is commonly used :

R(x) = x TV = p i=1 [∇x] i 2 , (II.18)
where [∇x] i is the finite difference approximation of the gradient of the continuous signal underlying x at x[i]. If x represents a 1D signal, [∇x] i is a scalar and the total variation closely resembles a Tikhonov regularizer with a Tikhonov matrix chosen a discrete gradient operator. However, one can get a good intuition of the key difference between the two by taking a bayesian standpoint. Indeed, one can interpret the total variation constraint as following from a probabilistic Laplace model on the gradient distribution of the underlying signal in the 1D case. Thereby, all else being equal, strong gradient values are more likely for this model than for the gaussian model associated with a Tikhonov regularization. Thus, sharp discontinuities can be allowed in the MAP solution while spurious variations associated with noise remain penalized.

As we will see in the next section, the discrete total variation can be seen as a particular instance of a larger class of non-smooth regularizers relying on deterministic low complexity priors on the signals. We complete these examples with ∞ norm based regularization. Indeed, in certain applications mostly related to digital communications, one might want to recover a flat signal or correlatively, compute and transmit a "democratic" or "anti-sparse" representation of a given signal for noise resilience and efficient amplifiers sizing (see the references and examples in [Elvira et al. 2015]). This is done by imposing a sufficiently low maximal amplitude bound on the signal to be computed, which can be achieved by choosing R(x) = x ∞ . The penalized regressions presented so far can be interpreted as MAP estimations yielding the so called bayesian methods ; in a statistical framework, the problem of choosing the regularization parameter is exchanged with that of modeling the uncertainties on the observed and underlying signals (see [START_REF] Wong | [END_REF] and the references therein). However, let underline the MAP interpretation might be misleading because -the MAP estimators do not follow in general the probability distribution used as a prior [Nikolova 2007] and -their might be other valid bayesian interpretations of the MAP estimators [Gribonval 2011].

Implicit regularizers We now give examples of regularization schemes that unlike the methods presented so far do not calculate the minimum of an explicitly defined functional. We first consider the Truncated Singular Value Decomposition (TSVD) method [START_REF] Hansen | [END_REF]]. We recall that the observation matrix M ∈ M pm (R) (Eq. II.8). The SVD of M is given by following factorization :

M = UΣV H (II.19)
where the matrices U ∈ M p (C) and V ∈ M m (C) are unitary and Σ is a (rectangular) diagonal matrices which diagonal elements are noted

σ 1 ≥ • • • ≥ σ min(p,m) ≥ 0.
The scalars σ i are called the singular values of M and somehow generalizes the notion of eigenvalues to non-square matrices. Thus, (II.20) where t = min(p, m). We note t + the index of the smallest non-zero singular value of M. The vector

M = t i=1 σ i U[:, i]V[:, i] H ,
x = t + i=1 U[:, i] H y σ i V[:, i] (II.21)
is the solution of the optimization problem

min x 1 2 x 2 2 s.t. x = argmin t 1 2 y -Mt 2 2 (II.22)
In words, x is the minimal energy minimizer of the functional J(x) = 1 2 y-Mx 2 2 ; x is therefore the unique solution when the problem is well-posed. If M represents a circular convolution, then V H is nothing but the matrix representation of the Discrete Fourier Transform (DFT) and U = V [START_REF] Gray | [END_REF]]. Thereby, U's columns are discretized complex exponentials at different frequencies. Yet, in deconvolution problems as previously illustrated by Fig.II.2, the matrix M is typically ill-conditioned and the lowest singular values are those associated with the highest frequencies. Therefore the high frequencies content of y, which is in general dominated by noise, is amplified in x. The TSVD then simply consists in dumping the terms in Eq.II.21 associated with the lowest singular values :

xTSVD = t SVD i=1 U[:, i] H y σ i V[:, i], (II.23)
for some user defined integer t SVD < t + . The conditions under which the TSVD yields similar results as the Tikhonov regularization in the case Γ = I p are studied in [START_REF] Hansen | [END_REF]]. One can note that this method makes no attempt to reconstructing the unobserved components of the image i.e. the projection of x on the null space of M.

Our second example is referred to as the early stopping in computational statistics and the machine learning literature (see [START_REF] Raskutti | Early stopping and non-parametric regression : an optimal data-dependent stopping rule[END_REF]] and the references therein). It relies on signals restoration using iterative schemes. Indeed, as shown in the previous examples, solving inverse problems is in general formulated as finding a stationary point of a certain cost function. While in several interesting cases there is no closed-form expression to compute the stationary points, in most cases, it can be done iteratively. We review some families of optimization methods in a forthcoming section. Iterative schemes generate a sequence (x k ) that converges in some sense toward a stationary point x of the (regularized) cost function, up to the computer's numerical precision. The restored signal is given by x = x kmax , where k max is set in situ, in general based on a numerical criterion of convergence. If the problem is not regularized, early iterates capture the low frequency content of the observed signal whereas as k increases, the iterates get corrupted by noise. Hence, the number of iterations itself behaves as a regularization parameter, which is the whole point of early stopping. This regularizing effect can be precisely characterized for the simple Landweber iterations scheme [Landweber 1951] which uses in the update formula :

x k+1 = x k -µM T (Mx k -y), (II.24)
where k ∈ N, µ ∈]0, 2/σ 2 1 [, σ 1 being the highest singular value of M. If x 0 = 0, then one can verify that the k th iterate is given by [Karl 2005]

x k = t + i=1 (1 -(1 -µ * σ 2 i ) k ) U[:, i] H y σ i V[:, i],
(II.25)

where t + has been previously defined as the index of the smallest non-zero singular value of M. Thus, the Landweber scheme effectively "re-weights" the un-regularized problem minimal energy solution given in Eq. II.21 with a set of varying weights identified as Loosely speaking, the features associated with the singular values smaller than σ are penalized while the features associated with the singular values larger than σ are kept. Thus, the number of iterations plays the role of the inverse of the regularization parameter in the previously described methods.

w i,k = 1 -(1 -µ * σ 2 i ) k , (II.
When the matrix M T M is invertible, the unregularized problem has a unique optimum, solution of the positive definite system

M T Mx = M T y.
(II.28)

Thus it can be calculated using the conjugate gradient (CG) algorithm [Hestenes & Stiefel 1952] which typically requires less iterations than the Landweber scheme. It has been shown that the k th iterate generated by the CG is the solution of the problem where

min x 1 2 y -Mx 2 2 s.t. x ∈ K k (M T M, M T y) (II.29)
K k (M T M, M T y) = span({M T y, (M T M)M T y, • • • , (M T M) k-1 M T y})
is the Krylov subspace generated by the matrix M T M and the vector M T y [Freund 1992]. Again, an early stopping in the CG method regularizes the inverse problem and the regularization weakens when the number of iterations (and consequently the dimension of the Krylov subspace) increases.

We finally note with no formal description that the early stopping is often used in the celebrated Richardson-Lucy to mimic regularization [Prato et al. 2012].

B Low complexity priors

In this section, we focus on the low complexity priors based regularizers. The general idea is that the estimated signal only needs a few number of parameters w.r.t. the ambient dimension to be fully specified, which can be used to recover well-posedness. In particular, we consider linear low complexity models [Baraniuk & Wakin 2009] which assume that the signal of interest can be explicitly written as a linear combination of few elements from a predefined or adaptive collection of atoms that will be referred to as a dictionary, following the terminology introduced in [Mallat & Zhang 1993]. This is described explicitly in the forthcoming sections.

B.1 Sparse models

B.1.1 Definition

A vector x ∈ R m is strictly sparse if most of its entries are equal to zero :

x 0 def = card(supp(x))
m.

(II.30)

x is said to be k-sparse if x 0 = k. Sparse signals appear in various applications. For instance in geophysics, the underground reflectivity can be modeled as a sparse signal (see [Claerbout & Muir 1973]). One can also think of acoustics and especially Room Impulse Response modeling [Sturm & Defrance 2010], and biological signal processing (for instance electromyography signals [START_REF] Liu | [END_REF]).

However, natural and especially physical signals are not strictly sparse in general. Yet, they might be compressible in the sense that most of their energy is captured by a few entries. This can be formalized in terms of the entries magnitudes decay's rate : assuming that σ m is a permutation in 1, m so that

|x[σ(1)]| ≥ • • • ≥ x[σ(m)]|, x is compressible or weakly sparse if ∀k ∈ 1, m , |x[σ(k)]| ≤ Rk -s , (II.31)
for some strictly positive reals R and s.

We note x (k) the best k-sparse approximation of x, for some k < m. x (k) is obtained by setting the m -k lowest amplitude entries of x to 0 so that the approximation error verifies

E x,k def = x -x (k) 2 2 = m i=k+1 x[σ(i)] 2 ≤ R 2 k -2s+1 2s -1 , (II.32)
if s > 1/2. Thus, for a given level of accuracy, the higher is s, the lower k can be chosen and therefore the more compressible is the signal. The definition above describes signals which are approximately sparse in their direct domain of representation, or more formally, in the canonical basis of R m . Yet, structured signals such as piece-wise smooth or oscillatory ones can be represented in alternative basis which better takes into account their specificity, yielding a faster decay of the sorted representation coefficients amplitudes. This is illustrated in the figures II.5 and II.6 for a gravitational waves signal (see [START_REF] Abbott | Observation of gravitational waves from a binary black hole merger[END_REF] for more details).

Therefore, more generally a signal is compressible if it can be written as

x = K i=1 a i d i , (II.33)

B.1.2 Promoting sparsity

Now we turn to the linear inverse problem of Eq.II.8 and we assume that the signal x is at least weakly sparse, w.r.t. some dictionary D. This prior can be naturally taken into account through the following constrained optimization problem :

min a 1 2 y -MDa 2 2 s.t. a 0 ≤ β, (II.35)
where β sets the minimal level of sparsity required. Let consider the set of admissible solutions, often referred to as the feasible set in numerical optimization literature :

B 0 R K (0, β) def = {a ∈ R K / a 0 ≤ β}. (II.36)
Assuming that β ≥ 1, it can be decomposed as follows :

B 0 R K (0, β) = β k=1 J∈P( 1,K )/card(J)=k span({d i } i∈J ), (II.37)
where for a set E, P(E) is the set of subsets of E. The structure of this set shows how difficult the resolution of the problem II.35 might be :

-B 0 R K (0, β) is non-convex because of the denumerable unions, while convexity is crucial in solving efficiently optimization problems ; the interested reader might have a look at [Boyd & Vandenberghe 2004] ; -as K increases, the number of subsets in B 0 R K (0, β) rapidly becomes too large for a brute force global resolution of the problem II.35 to be tractable.

Thus, the 1 norm is generally used as a convex relaxation of the 0 pseudonorm, which have been justified in different ways. The most natural argument is based on the simple remark that the unitary

1 ball B 1 R K (0, 1) is the convex hull of the intersection of the unitary 0 ball B 0 R K (0, 1) with the unitary ∞ ball B ∞ R K (0, 1).
In [START_REF][END_REF], Linear Programming (LP) arguments are used to support the fact that 1 minimization promotes sparsity. Thereafter, outstanding results are established showing that if the underlying signal is sufficiently sparse (or compressible), 1 norm minimization can yield an exact or stable recovery [START_REF] Elad | [END_REF], Donoho 2006, Candès et al. 2006]. More recently, a generic and formal framework has been introduced in [START_REF] Chandrasekaran | [END_REF] for translating simplicity priors into convex penalties ; in particular, the 1 norm emerges in this framework as the atomic norm associated with the set of 1-sparse unitary vectors in the sparse signal's space.

Hence, we consider the optimization problem

min a 1 2 y -MDa 2 2 s.t. a 1 ≤ β, (II.38)
or equivalently in its augmented lagrangian form

min a 1 2 y -MDa 2 2 + λ a 1 (II.39)
for some well-chosen λ [START_REF] Ciak | [END_REF]]. These problems can be solved efficiently.

In certain applications such as DNA microarray analysis [Schena et al. 1995], the signal of interest is known to have groups of entries which are correlatively significant or negligible. The paradigms of group or block sparsity have been introduced for taking into account this particular prior [Huang et al. 2010, Peyré & Fadili 2011]. Specifically, one needs to define a partition

{G 1 , • • • , G q } of 1, K which we recall requires that q i=1 G i = 1, K and G i G j for i = j.
The ell 1 norm is then replaced by the 1,2 norm defined as

a 1,2 def = q i=1 a (G i ) 2 , (II.40)
where a (G i ) is the vector obtained by extracting a's entries indexed by G i .

Remark : note that the earliest sparse recovery attempts consisted in solving approximately the 0 constrained problem II.35 via greedy methods, namely the celebrated Matching Pursuit [Mallat & Zhang 1993] and the Orthogonal Matching Pursuit [Pati et al. 1993]. Several methods have been proposed since then for 0 minimization, like for instance the Iterative Hard Thresholding (IHT) (see for example [START_REF] Starck | [END_REF], Starck et al. 2005]) ; the IHT is thoroughly analyzed in [Blumensath & Davies 2008] and [Blumensath & Davies 2009]. Let finally note a more recent trend which consists in approximating the 0 pseudo-norm with the so called p norms, with 0 < p < 1 [START_REF] Zheng | Does p -minimization outperform 1 -minimization ?[END_REF].

B.1.3 Analysis formulation of the sparsity prior

According to the model II.33, the coefficients computed by solving II.39 enable one to (approximately) synthesize the underlying signal x using the dictionary D atoms. Hence the terminology of synthesis formulation used when D = I p . Alternatively, one can adopt the analysis formulation which has the following form :

min x 1 2 y -Mx 2 2 + λ D T A x 1 , (II.41)
for some dictionary D T A . In this case, one looks for a signal compatible with the data and strongly correlated with only a few atoms in D A . Although very similar, the analysis and synthesis formulation are only equivalent if the dictionary involved is square and invertible [Elad et al. 2007]. Precisely, if D T A = D -1 , then a vector x is the solution of Problem II.41 if and only if D T A x is the solution of Problem II.39. This is crucial from applications standpoint because of the importance of overcomplete dictionaries in sparse recovery as we will shortly detail. Let first mention that the total variation based regularization falls into this framework in the 1D case with D A chosen as the transpose of the matrix representing of the finite difference approximation of the gradient. Therefore, this penalty is particularly suitable for (approximately) piece-wise constant signals or equivalently signals that have a sparse gradient. The total variation in 2D can also be expressed in the analysis framework, but involving an 1,2 norm (see Equation II.40).

In this example, the analysis formulation is weaker that the synthesis one in a certain sense. Indeed, when the dictionary D A has fewer atoms than the underlying signal's space dimension, the analysis enables one to promote some regularity on the signal without completely specifying its structure. Precisely, the signal is only imposed to be orthogonal to certain columns of D A whereas in synthesis, the signal is explicitly described as a linear combination of some of the dictionary's atoms.

B.1.4 The dictionary

The dictionary choice is in first instance driven by the signal's structural specificity. Indeed, the signal can only be sparsely encoded in the dictionary if some of its atoms closely resemble the signal main features. Let give a few classical examples :

-the DCT dictionary is suitable for oscillatory signals as already illustrated ;

-standard wavelet bases are optimal for encoding signals with isolated singularities [Mallat 1999] ;

-in presence of "extended" singularities with geometrical regularity such as edges, the curvelet dictionary are more adapted [Candès & Donoho 2004].

Because of the tremendous morphological diversity of natural signals and the success of sparse recovery in various applications (see [Zhang et al. 2015] and the references therein), the task of designing increasingly adapted transforms and dictionaries has considerably polarized the recent signal processing literature. This is well illustrated by [Jacques et al. 2011] and [Duval 016]. One can note that a particular focus has been put on the so-called multiscale representations. Indeed, characteristic features can appear at multiple scales, which in general are not known beforehand. Similarly, translation invariance is an important dictionary feature [Coifman & Donoho 1995]. Besides images generally contain structurally different features, so that each "structural" component of the image can be sparsely represented in a different dictionary. In this case, a better dictionary, from the sparse recovery standpoint, is obtained by concatenation. This is the basic rational behind the Morphological Component Analysis [START_REF] Starck | [END_REF]] and its different extensions. More generally, geometrical invariances requirement and transforms combination yield overcomplete dictionaries, which have been proven very powerful in sparse signal recovery [START_REF] Starck | [END_REF]. We now recall two theoretical results which are important to bear in mind when building or choosing a dictionary. The first one concerns the uniqueness of sparse representations in overcomplete dictionaries. We consider two orthonormal bases

D 1 = [d 11 , • • • , d 1p ] and D 2 = [d 21 , • • • , d 2p ]
, and we consider a vector x in R m we note α the sparsest vector verifying

x = [D 1 , D 2 ]α.
(II.42)

It is shown in [START_REF] Elad | [END_REF]] that if

α 0 ≤ 1 2 (1 + 1 µ(D 1 , D 2 )
), (II.43) then α is unique and can be calculated via an

1 minimization, µ(D 1 , D 2 ) = max 1≤i,j≤p
{d T 2i d 1j } being the so-called mutual coherence of the basis D 1 and D 2 . This suggests that in general, having strongly correlated atoms withing the dictionary is susceptible to undermine the ( 1 norm based) sparse recovery. The second result gives an indication on the dictionary choice, with respect to the observation operator in Problem II.39. Precisely, we assume, among other omitted hypothesis, that there exists an m × m -p matrix U so that M = [M T , U] T is an orthogonal matrix. Then, the smaller is µ( MT , D), the less observation are needed for achieving a given level of accuracy in the sparse recovery, in probability [Candès & Romberg 2007].

In words, the more D atoms have spread representations w.r.t. MT , the more the recovery is robust to under-determination.

Dictionary learning So far, we have mentioned dictionaries derived analytically and which are in general optimal for representing certain formal classes of signals. However, in certain applications, one disposes of data samples that are structurally similar to the signal. Hence the natural idea of learning the dictionary directly from the data. Precisely, given a training set x 1 , • • • , x q , the core idea is to simultaneously calculate a dictionary D L and a set of sparse vectors a 1 , • • • , a q so that x i ≈ D L a i , (II.44)

∀i ∈ 1, q . Dictionary learning has been successfully applied in various applications like for instance in astronomy [Beckouche et al. 2013b], for Blind Source Separation [Bobin et al. 2007], neuroscience [Varoquaux et al. 2011] and more generally natural images restoration [START_REF] Elad | [END_REF], [START_REF] Mairal | [END_REF]. It seems worth noting that while most dictionary learning methods focus on the synthesis formulation of sparsity, a few work exists that are based on the analysis model (see [Rubinstein et al. 2013] and the references therein).

In both cases, the dictionaries learned are in general unstructured so that they have to be fully stored and explicitly used in matrix-vector products. This is an important consideration from a practical point of view since non-adaptive dictionaries such Wavelets, DCT or Curvelets dictionaries are associated with fast implicit implementations, which makes them computationally efficient.

B.2 Low rankness

In this section, the signal of interest is represented by an m × q matrix X = [x 1 , • • • , x q ] which is convenient in various applications such as video or multichannel signal processing. The observation model is adapted as follows :

Y = M(X) + B, (II.45)
where M is a linear operator and Y ∈ M pq (R).

B.2.1 Definition

The (algebraic) rank is defined as follows :

rank(X) def = dim(span((x i ) 1≤i≤q )) = rank(X T ).
(II.46)

X will be said to have a low rank if rank(X) min(m, q). (II.47)

To get an intuition of the practical interest of this quantity, let assume that X is an hyperspectral image, with the convention that the i th column x i corresponds to the i th band image. Each line of X can be modeled as a linear combination of deterministic spectra related to the materials present in the scene [Chang 2003]. Thus, if the scene contains exactly s different materials, then the rank of X is at most equal to s ; this is particularly interesting when only a few materials are preponderant in the scene. As illustrated by this example, in general, constraining the rank allows one to account for the redundancies present in the data to be restored. We note a(X) = (a i ) 1≤i≤t the singular values of X, with t = min(m, q). Then we have rank(X) = a(X) 0 .

(II.48)

Recalling that X can be written in the form

X = t i=1 a i U i , (II.49)
where the U i s are rank one matrices, the low rankness clearly appears as an extension of the notion of compressibility to matrices. Precisely, the rank quantifies the complexity of a linear model that would perfectly fit the data stored in the considered matrix.

B.2.2 Constraining the rank

As seen in the previous section for the 0 norm, the direct problem

min X Y -M(X) 2 F s.t. rank(X) ≤ β (II.50)
is combinatorial and NP-hard [Fazel 2002]. Thus, various heuristics have been proposed for promoting low rankness.

Trace heuristic The trace function has been successfully used as a surrogate for the rank [Pare 2000] when the estimated matrix is positive semidefinite matrix (PSD). This yields a convex tractable problem. Observing that in this case Trace(X) = a(X) 1 (see Eq.), this heuristic can be understood as an 1 based sparsity promoting approach.

Log-det heuristic In a similar vein, a particularly fruitful idea has been to use the function Y → log(det(Y + δI m )) as a smooth surrogate of the rank function on the cone of m × m PSD matrices [START_REF] Fazel | [END_REF]]. While being concave, it yields a reweighted-1 like minimization on the singular values [Candès et al. 2008]. Thus it can be considered as a refinement of the trace heuristic. In the same work, the authors derive an extension to arbitrary non square matrices.

Nuclear norm heuristic The most popular popular rank surrogate function is certainly the nuclear norm defined as the sum of the singular values of a given matrix [START_REF] Fazel | [END_REF]] :

X = a(X) 1 . (II.51)
This extends the trace heuristic to arbitrary non square matrices. The nuclear norm can be considered the most suitable convex function for promoting low rankness in more than one respect. Similarly to the 1 norm, the nuclear norm has been shown to be the convex envelop of the rank function on the set matrices {Y/ S(Y) ∞ ≤ 1} [START_REF] Fazel | [END_REF], where S(Y) is the vector of singular values of Y and therefore, S(Y) ∞ is Y's spectral norm. A more recent point of view, previously invoked for the 1 norm, derives the nuclear norm as the atomic norm associated with the set of rank one matrices with unitary spectral norm [START_REF] Chandrasekaran | [END_REF]. Most importantly, various low rank matrices recovery guarantees using nuclear norm have been established, with a striking analogy to 1 based sparse recovery results (see [Recht et al. 2010, Candès & Plan 2010, Candès et al. 2011, Candès & Recht 2009]).

Low rank minimization has been successfully applied in various task such as matrix completion, for recommendation systems or missing data problems, background detection in videos and hyperspectral images restoration.

Matrix factorization

In certain applications, the rank of the matrix under estimation that we denote r is known. While it is not straightforward to translate this information into a suitable value of the regularization parameter in a nuclear norm minimization, one can directly impose the rank through a matrix factorization. Indeed, the following equivalence holds :

rank(X) ≤ r ⇔ ∃(S, A) ∈ M mr (R) × M rq (R)/X = SA (II.52)
which suggests the alternative formulation of the inverse problem

min (S,A)∈Mmr(R)×Mrq(R) Y -M(SA) 2 F . (II.53)
This problem is globally non convex. However it is bilinear in S and A and hence partially convex for each of these matrices. Moreover this approach enables one to incorporate more priors by adding constraints on S and A. This yields important frameworks such as Non-negative Matrix Factorization (NMF) [Wang & Zhang 2013], or Sparse NMF [START_REF] Kim | [END_REF].

Matrix factorization and nuclear norm minimization have complementary aspects and depending on the prior knowledge available, one might be more suitable than the other. However, an attempt of using the two in an unified framework can be found in [Cabral et al. 2013].

Remark : in certain applications such as system reduction or identification, the low rank matrices of interest are known to have a specific structures which can be usefully integrated in the estimation process. This yields alternative rank minimization methods and algorithms [Markovsky 2008]. In Chapter II, we presented how signals structural properties can be used for regularizing inverse problems, with a focus on linear low complexity data models. This chapter goes beyond the framework of dictionaries based representations and transforms and considers more general non linear data parametrizations ; the main idea being to uncover and exploit the non linear geometry underlying a data set, especially when this geometry is simple, in a sense to be defined.

Chapitre III

Beyond linear data models

In Section A.2, we recall important geometric notions with this respect. We refer the interested reader to [Flaherty & do Carmo 2013] for a more detailed presentation. Then we present an overview of Manifold Learning (ML) main techniques. Section B is an introduction to Optimal Transport (OT). The emphasis will be put on certain geometrical properties of OT that are particularly interesting from a data analysis standpoint.

A Manifold data model

A.1 Reminder on differential geometry

A.1.1 General definitions

A d-dimensional manifold can be though of as a collection of points M which resembles R d in each points neighborhood. For instance an 2 sphere in R 3 is a 2-dimensional manifold. We precise the definition using the notions of charts and atlases. A bijective map φ between a subset X of M and an open subset of R d is called a d-dimensional chart of M denoted (X , φ). For each point x, the vector φ(x) constitutes its coordinates in R d ; φ -1 is referred as a parametrization of M around x. Thus, M can be viewed as an union of suitably chosen coordinates patches or charts.

It follows naturally that each point of M shall belong to at least one chart's domain and that two overlapping charts shall give consistent information, which yields the notion of atlas.

A smooth atlas of M in R d is a collection of d-dimensional charts (X i , φ i ) i∈Ω for some indexation domain Ω verifying 1. i∈Ω X i = M 2. for (i, j) ∈ Ω 2 , i = j/ X i X j = ∅, then the sets φ i (X i X j ) and φ j (X i X j ) are open and the map φ i • φ -1 j : R d → R d is smooth on φ j (X i X j ).
If M is a set of geographical locations on the Earth surface, then the notion of atlas recovers its common assertion.

The maximal atlas A + generated by an atlas A is the atlas defined as the set of the all the charts (X , φ) such that A (X , φ) is also an atlas of M. This can be viewed as the minimal structure required for extending differential calculus from linear subspaces to arbitrary sets. Indeed, the study of a functional f : M → R can be replaced on the domain of a particular chart (X , φ) in A + with the study of the vector space functional

f : R d → R, x → f (φ -1 (x)), (III.1)
on the set φ(X ). Thus, a maximum atlas is also referred as a differentiable structure on the M. From this, a manifold can be formally defined as couple (M, A + ) of a set of points and its differentiable structure (when it exists). By a slight abuse of notation, in general we will only use M for referring to the manifold. There are numerous interesting examples of manifolds. At this point we only mention the particular case of vector spaces which, thereby, are often referred as linear manifolds, the subject of this section being the non linear manifolds. In the following, we focus on connected manifolds i.e. manifolds which points can not be represented as the union of two or more disjoint nonempty open subsets.

A.1.2 Manifold's first order geometry

A curve on a manifold M is defined as a function γ : R → M. Assuming that M is a submanifold of a vector space E, a tangent vector at a point p of

M is any vector v in E such that v = lim t→0 γ(t) -γ(0) t def = γ (0), (III.2)
for some smooth curve γ on M which verifies γ(0) = p. The tangent space to M at p denoted T p M is the collection of all the tangent vectors at p. If M is a d-dimensional manifold, one can verify that T p M is a vector space of dimension d. We note that d is often referred in applications as the intrinsic dimension of M by contrast to the dimension of the ambient space E. As we will see in Section A.2, this notion is important from a practical point of view, especially when d dim(E). T p M can be endowed with an inner product ., . M,p which induces the norm

x M,p = x, x M,p . (III.3)
The manifold M is Riemannian if ., . M,p varies smoothly with p ; in this case, the inner product ., . p is referred as the Riemannian metric.

Riemannian manifolds and more generally, the Riemannian geometry constitute the mathematical framework of the General Relativity ; precisely the theory relates the Riemannian metric ., . M,p also referred as the metric tensor at every space-time point to the mass/energy distribution through the celebrated Einstein's equations [Besse 2007]. In this physical context, the smoothness requirement is particularly adapted.

In the Riemannian framework, the length of a smooth curve γ

: [t 1 , t 2 ] → M is given by L(γ) = t 2 t 1 γ (t) M,γ(t) dt. (III.4)
This gives a mean of measuring distances on the manifold. Indeed, one can introduce the set of smooth curves joining two points p 1 and p 2 in M :

Γ(M, p 1 , p 2 ) = {γ : [0, 1] → M/ γ(0) = p 1 , γ(1) = p 2 and γ smooth on [0, 1]}.
(III.5) Then the length of the shortest path between p 1 and p 2 can be calculated as

L M (p 1 , p 2 ) = argmin γ∈Γ(M,p 1 ,p 2 ) L(γ).
(III.6)

The functional L M : M 2 → R known as the Riemannian distance actually defines a metric on M, which makes a connected Riemannian manifold a metric space.

If we consider a d-dimensional chart (X , φ) of M, φ is a local isometry of (X , L M ) in R d , 2 if ∀ (p 1 , p 2 ) ∈ X 2 L M (p 1 , p 2 ) = φ(p 1 ) -φ(p 2 ) 2 . (III.7)
Such property is particularly suitable in dimension reduction tasks when the data are sampled from a Riemannian manifold. Furthermore, if L M is a euclidean norm, then the manifold is qualified as intrinsically euclidean.

A.1.3 Geodesics on Riemannian manifolds

Certain curves on M have the property that they induce shortest paths locally. It is the case of the great circles if M is a sphere in R 3 . These curves are called geodesics and we recall thereafter some related properties [Lee 2006]. Formally, a curve γ on a Riemannian submanifold M of a vector space is a geodesic if the acceleration vector is permanently normal to the tangent plane :

γ (t) ∈ T γ(t) M ⊥ (III.8)
for any t in the domain of γ. This implies in particular a constant velocity along the curve :

γ (t) M,γ(t) = c. (III.9)
Recalling Eq.III.4, this yields an interesting property of geodesics, which is, they are parametrized with an affine function of the arclength. As a simple example, the geodesics in a vector space are the constant speed straight curves.

A geodesic is minimal if its length is smaller than the lengths of all the smooth curves on the manifold joining its endpoints. The length of the minimal geodesic is nothing but the Riemannian distance between its endpoints, hence the terminology of geodesic distance sometimes found in the literature. Any sufficiently small segment of a geodesic is minimal. An interesting property from a practical point of view is that if a geodesic is minimal on an interval I, then its restriction to any subinterval of I is also minimal. However, the shortest path between two points of a connected Riemannian manifold need not to be a geodesic, even up to a change of parametrization. Geodesically complete Riemannian manifolds are precisely those in which any two points can be connected through a geodesic. In particular, in this case minimal geodesics exist for any two points.

A nice connection between the metric space (M, L M )'s topology and M's geodesics is given the Hopf-Rinow theorem [Hopf & Rinow 1931] : among other equivalences, it states that (M, L M ) is (topologically) complete if and only if M is geodesically complete.

A.2 Manifold learning techniques

The computational applications of manifold data models are two-fold. On the one hand, there are situations where one needs to estimate a signal that belongs to a known manifold. One can find some examples related to matrices manifolds in [START_REF] Absil | [END_REF]]. An application to natural signals and images restoration can be found in [START_REF] Peyré | Manifold models for signals and images[END_REF]]. This first framework boils down to solving optimization problems on non linear manifolds, which constitutes a research area on its own. Among other recent contributions, we can cite [Bacák et al. 2016] and some references therein. Conversely, one might have a data set (approximately) sampled from an unknown manifold and instead try to uncover the underlying manifold geometry. This can serve purposes of visualization, classification or data synthesis. We term Manifold Learning (ML) the corpus of techniques aiming at characterizing a manifold from a set of samples lying on or close to the manifold. This has been a particularly active field in the Machine Learning community over the past decade. We review thereafter some of the most influential contributions.

Unless mentioned otherwise, the following notations holds in this section :

-(x i ) 1≤i≤N is the learning data set ;

-the learning set is sampled from a submanifold M of R m ; -(y i ) 1≤i≤N is the low dimensional representation of the training set ;

d might refer to the intrinsic dimension of M or the dimensionality of (y i ) 1≤i≤N .

A.2.1 Multidimensional scaling

The multidimensional scaling (MDS) refers to a set of techniques that aims at determining the underlying attributes and the spatial organization of a set of objects based on measurements of pairwise similarities (or dissimilarities) between the objects. We limit our presentation to the so-called classical MDS. This method was originally introduced in Psychometrics for uncovering stimuli perceptual dimensions and their scales [Torgerson 1952]. Given a set of N stimuli, one is given a matrix D = (d ij ) 1≤i,j≤N , where the scalar d ij is a comparative distance between the i th and the j th stimuli, which we assume here to be positive. However, this matrix does not need to derive from a metric function. The classical MDS computes a set of vectors

(y i ) 1≤i≤N in R d so that y i -y j 2 ≈ d ij , ∀ (i, j) ∈ 1, N 2 by solving the optimization problem min y 1 ,••• ,y N (i,j)∈ 1,N 2 (d ij -y i -y j 2 ) 2 . (III.10)
We note that this Problem has a known closed-form solution based on an eigen-system decomposition of D ; this is presented in more details in Chapter VI. The dimension r is a parameter and the coordinates have hierarchical importance, which follows directly from the solution's structure. This allows one to visualize the object space at different scale (when d > 3). The MDS has come to be a popular tool for exploratory data analysis. A straightforward application is the problem of recovering city coordinates from pairwise distances. Although it is not directly meant for manifold-valued data analysis, MDS constitutes an important step in the foundational ML technique Isomap, which is presented in the following section. We refer to [Borg & Groenen 2005] for a thorough presentation of MDS.

A.2.2 Isomap

The Isomap is a non-linear dimension reduction method introduced in The dimension parameter r introduced in the previous section shall be ideally chosen as the intrinsic dimension of M. It corresponds to the number of degrees of freedom of the system that generated the data. Isomap is guaranteed to recover the true underlying geometrical structure as the sampling density tends to the infinity if the manifold is intrinsically euclidean. Besides, it can be successfully applied for several images manifolds of interest in computer vision [Donoho & Grimes 2005]. However, the Riemannian distances estimation step might be computationally demanding, since it involves finding shortest paths on a potentially large weighted graph. This shortcoming is avoided in the next approach.

A.2.3 Local Linear Embedding

We consider the same setting as in the previous section. The Local Linear Embedding (LLE) [Roweis & Saul 2000] lies on natural geometrical intuitions. Firstly, if the samples are sufficiently dense then each sample's nearest neighbors approximately lie in its tangent space (in fact an affine space parallel to this tangent space). Hence, each sample can be expressed as a linear combination of its nearest neighbors. Hence, the first step of LLE is to compute weights matrix W = (w ij ) 1≤i,j≤N by solving the following problem :

min W N i=1 x i - j∈nn k (i) w i,j x j 2 2 , s.t. j∈nn k (i) w i,j = 1 ∀i ∈ 1, N , (III.11)
where nn k (i) is the set of indexes of the k nearest neighbors of the sample

x i . Besides w ij = 0 if j / ∈ nn k (i).
These weights are invariant under isometric affine transformations and rescaling of the manifold. Thereby, they locally capture its intrinsic geometry, and should remain optimal in a low dimensional representation of the manifold that "preserves" this geometry. It follows that the low dimension coordinates can be estimated by solving

min y 1 ,••• ,y N N i=1 y i - j∈nn k (i) w i,j y j 2 2 , s.t. j∈nn k (i) w i,j = 1 ∀i ∈ 1, N . (III.12)
Interestingly, this method derives a global characterization of the manifold from local information, which is made possible by the neighborhoods overlaps. The same principle will appear in the next methods.

A.2.4 Laplacian Eigenmaps

As for the LLE, the Laplacian Eigenmaps methods builds over the idea of preserving local neighborhood information. We note nn(i) the set of indexes of the "relevant" neighbors of the i th sample. The low dimensional coordinates are calculated by solving the problem

y 1 ,••• ,y N N i=1 j∈nn(i) w ij x i -x j 2 2 , (III.13)
where

w ij = exp(- x i -x j 2 2 t ), (III.14)
with some t > 0. Thus, neighbor samples are strongly constrained to remain close in the new coordinate system. As previously, let introduce the weights matrix W = (w ij ) 1≤i,j≤N , with the convention

w ij = 0 if j / ∈ nn(i). It is
shown in [Belkin & Niyogi 2003] where D is a diagonal matrix defined as D[i, i] = j∈nn(i) w ij , for i ∈ 1, N and L = D-W. L is the laplacian matrix of a weighted graph G with N nodes representing the data sample, the weight of the vertex between the i th and the j th node being equal to w ij [START_REF][END_REF]]. Interestingly, it is shown that for t sufficiently small, L can be considered a discrete approximation of the Laplace-Beltrami operator on M over the graph G. Besides, the eigenfunctions of the Laplace-Beltrami operator on M are solution of the problem

min f M ∇f (x) 2 2 dµ(x), s.t. M f (x) 2 dµ(x) = 1 (III.16)
for some probability measure dµ on M. This shows that the Laplace-Beltrami operator's eigenfunctions tend to realize mappings of the manifold that moderately affect the local samples organization. This supports the geometrical intuition which underpins the Laplacian Eigenmaps method.

A.2.5 Hessian Eigenmaps

The Hessian-based Local Linear Embedding (HLLE) methods represents an improvement over each of the pre-cited ML techniques in terms of preservation of the manifold's (local) geometry [Donoho & Grimes 2003]. At the heart of this method, a quadratic form termed the H-Functional which is defined as

H(f ) = M H f (x) 2 F dµ(x), (III.17)
where f is a smooth functional on the manifold, H f (x) is the hessian matrix of f at x and dµ is a probability measure on M. The hessian matrix is defined in an orthogonal coordinates system in the tangent space T x M relatively to the dot product inherited from the ambient space. This quadratic form is well defined in the sense that it is independent of the orthogonal coordinates system chosen in a given tangent plane. It quantifies how curvy is a functional on the manifold, in average. To gain a better grasp of the intuition behind this choice, let consider the family of functionals (f i ) 0≤i≤m defined as

f i : R m → R, x → x[i], (III.18) ∀i ∈ 1, m and f 0 : R m → R, x → 1. (III.19)
These functions span the set of affine and hence, flat functions on R m . Similarly, for a non-linear manifold M, one can expect the set of flat functions on the manifold to be related to a canonical coordinate system on the manifold.

Precisely, if M is a Riemannian d -dimensional submanifold of R m which is locally isometric to an open connected subset of R d , then the null space of H is a d + 1-dimensional vector space consisting of the constant functionals on the manifold and a d-dimensional space spanned by the isometric coordinates. Thus, the HLLE derives the low dimensional representation by estimating the null space of H based on the learning set as follows :

-for each i ∈ 1, N , one computes a d(d + 1)/2 × k matrix H (i) from local tangent coordinates so that for a functional f defined over the learning set which is the discretization of a smooth functional f defined over M, if we note f (i) its restriction to a vicinity of size k of the i th sample, the vector H (i) f (i) constitutes an estimate of the entries of H f (x i ) ;

-from the matrices H (i) , one builds a N × N symmetric positive matrix H so that H(f ) ≈ f T Hf ;

-ones computes d + 1 orthonormalized eigenvectors corresponding to the d + 1 smallest eigenvalues of H ; we note these vectors u 0 , • • • , u d , where there are sorted in increasing order of the associated eigenvalues ; u 0 is associated with the eigenvalue 0 and corresponds to constant functionals ; the low dimension coordinates of the i th sample are given by the vector (u j [i]) 1≤j≤d .

The requirement for local isometry brings this method closer to the LLE in which it is implicit. Besides, the H-Functional somehow resembles the quadratic form minimized in Problem III.16, which can be rewritten in terms of the Laplace Beltrami operator. However, the hessian is somehow more suitable than a laplacian for linear functions recovery [Donoho & Grimes 2003]. Furthermore, the HLLE can be applied to a wider range of manifolds than Isomap which requires the manifold to be globally isometric to a convex subset of R r .

Remarks : in this thesis, we are interested in the problem of synthesizing new samples from a manifold (namely a PSFs manifold) given a learning sample set and according to a prescribed parametrization of the manifold. However most of the ML techniques previously presented do not provide a natural way of doing so, except for the LLE. For instance, in [START_REF] Tenenbaum | Joshua B Tenenbaumet al. Mapping a manifold of perceptual observations[END_REF]], generalized radial basis functions are used for computing a mapping between the low dimensional coordinates and the manifold. However, the recently introduced Geometric Multi-Resolution Analysis [Allard et al. 2012] seems particularly suitable in this regard. On the other hand, the reverse problem of determining the low dimension coordinates of an out-of-sample point has been treated in several works (see for instance [Bengio et al. 2003]).

For most manifold learning techniques, including those presented here, the manifold intrinsic dimension is a free parameter. Although in general, it can be empirically (and efficiently) tuned from the learning procedure itself, several methods have been proposed for estimating the intrinsic dimension of an embedded manifold, given a learning sample (see [Little 2011] and the references therein).

B Optimal transport

The photograph of the figure III.2 was taken by the author at the french "Japan Expo" and shows the display table of a mangas seller. Speculatively, we suppose that at the end of the exposition, the unsold books could fit in a single box. The optimal transport problem arises from two trivially observable facts in this situation :

-the seller can store the books in the box following different strategies ; for instance he can distribute each book pile in the box as evenly as possible, or he can try to preserve the book piles as much as possible ; -some strategies involve more efforts in terms of mass displacement than others. This naturally yields the question of finding the minimal effort strategy of transportation. The first known formal study of this question traces back to the late 18 th century ( [Monge 1781]). Since then, the Optimal Transport has come to be a remarkably fruitful topic in applied mathematics. In this section, we review the Optimal Transport's theory basic formalism and some important properties of the optimal solutions.

B.1 Problem formulations

B.1.1 Framework and basic definitions

Let consider two probability measures µ and ν defined respectively on some measure spaces X and Y . We shall define a measurable cost function c : X × Y → R + which quantifies the effort of transporting one unit of mass from a point in X to a point in Y . For a measurable map T :

X → Y and a subset B of Y , T -1 (B) def = {x ∈ X/ T (x) ∈ B}. (III.20)
We say that µ is the push-forward of ν by T if

for any measurable set B ⊂ Y, ν(B) = µ(T -1 (B)), (

B.1.2 Monge problem

The Monge Optimal Transport problem is the following :

min T X c(x, T (x))dµ(x), s.t. ν = T #µ (III.23)
where c(x, T (x))dµ(x) is by construction the cost of moving the elementary mass dµx from x to T (x). In several applications, X = Y and c is chosen as a metric on X. c is then termed the ground metric, in reference to the logistic interpretation of the Optimal Transport.

The Monge's problem imposes a transport with no mass split. In other words, the elementary mass dµx has to be mapped to a single location in Y . It appears intuitively that depending on µ and vu such a mapping might not be possible (for instance if X = Y = R and if one of the two probability measures has a dense support while the other has a discrete support in R). We will shortly introduce a relaxation in this regard due to Kantorovitch [START_REF] Kantorovitch | ON THE TRANSLOCATION OF MASSES[END_REF]].

B.1.3 Kantorovich problem

A transference or transportation plan between µ and ν is a probability measure Π on X × Y that verifies is a transportation plan between µ and ν so that P(µ, ν) is always non-empty. We refer the reader to [Villani 2003] for detailed discussions on the existence of solutions to the Monge and Kantorovich problems. However, we note that if X = Y = R m , for some m and if the cost is in the form c(x, y) = xy p p , (III.27)

for some p > 1, the Kantorovich problem in general admits a unique solution.

B.1.4 Brenier-Benamou problem

We assume that X = Y = R m . Let note T * a solution to the Monge problem whenever it exists. Provided that T * is smooth, µ and ν obey the relation µ

(x) = ν(T (x))| det(J T * (x))|, (III.28)
where J T * (x) is the jacobian matrix of T * at x [Ambrosio & Savaré 2007]. Let (T t ) 0≤t≤1 be a family of maps defined as

T t (x) = tx + (1 -t)T * (x), ∀(t, x) ∈ 0, 1 × X. (III.29) If c(x, y) = x -y 2 2 , then the function χ * defined on 0, 1 × X as χ * (t, x) = ν(T t (x))| det(J Tt (x))| (III.30) is solution of the problem min χ,v 0,1 ×X χ(t, x) v(t, x) 2 2 dtdx s.t. (III.31) ∂ t χ + div x (χv) = 0 (III.32) χ(0, .) = µ, χ(1, .) = ν, (III.33)
for a vector field v : 0, 1 × X → X. This is the celebrated Benamou-Brenier formulation of the Optimal Transport problem [Benamou & Brenier 2000]. The map t → χ * (t, .) can be viewed as a curve in the space of probability measures on X which joins µ to ν. This particular curve will be characterized more precisely in the next section. Let just say for now that it realizes the displacement of elementary mass particles of the distribution µ toward the distribution ν. Thus, v shall be interpreted as a velocity field and Eq.III.32 ensures mass preservation through the displacement. In the same spirit, χ(t, x) v(t, x) 2 2 has to be identified as a kinetic energy density in this formulation.

B.2 Geometry of optimal transport

B.2.1 The transport map structure

We consider a taxi firm which has three cars available located at x i , i ∈ {1, 2, 3} and needs to pick three potential customers up as fast as possible. The customers are waiting at y i , i ∈ {1, 2, 3} as depicted in Fig. III.4. Assuming homogeneous traffic conditions, this situation can be modeled as a particular instance of the Monge problem, opportunely using the Manhattan distance as the cost function.

The optimal assignment is easily picked up by eye. It also appears that

∀(i, j) ∈ {1, 2, 3} 2 , (x i -x j ) T (y σ (i) -y σ (j) ) ≤ 0 (III.34)
where σ (i) is the index that corresponds to the customer assigned to the i th cab. This seemingly coincidental fact illustrates a general property of the optimal transport solution, which turns out to be a monotonic rearrangement of one distribution onto the other under certain conditions [Villani 2003].

Precisely, using the same notations as in the previous section, with X = Y = R m and c(x, y) = xy 2 2 , if we further assume that µ and ν do not give zero probability to small sets (in terms of Hausdorff dimension), then there exists a convex function φ so that the optimal transport map is given by T = ∇φ.

(III.35)

Figure III.4 -Optimal (cab) assignment.
This result is often referred as the Optimal Transportation theorem and has been generalized to non quadratic costs. The previously mentioned monotonicity follows from the convexity of φ.

The existence of φ suggests potential applications of this framework to astronomical imaging. For instance, in strong gravitational lensing analysis, φ could be interpreted as the gravitational potential induced by a massive unknown structure.

We now turn to the geometry related to the minimal transportation cost.

B.2.2 Monge-Kantorovich distances

The Transport minimal cost in Monge and kantorovich problems measures the similarity between two distributions. Indeed, the more two distributions look alike, the less effort is needed for transforming one into the other, and conversely. This section precises this intuition.

Somehow the structure required for this purpose is weaker than what was previously needed in for characterizing the optimal transport map. In particular, X is only assumed to be a separable complete metric space, endowed with a metric D ; also Y = X. The cost function takes the following form : c(x, y) = D(x, y) q , (III.36)

for some q ≥ 1. Following the notations in [Villani 2003],

-T q (µ, ν) is the minimum cost associated with the Kantorovich problem ;

-P q (X) is the set of probability measures with finite moments of order q.

Then the functional W q : P q (X) × P q (X), (µ, ν) → T q (µ, ν)

1 q (III.37) defines a metric on P q (X). This metric is known under different names among which the Monge-Kantorovich distance of order q or the q th Wasserstein distance. In the particular case q = 2, it is termed the quadratic Wasserstein distance.

We assume that X = R m . P 2 (R m ) is a submanifold of the vector space of finite signed measures on R m . For a given µ ∈ P 2 (R m ), let ρ 1 and ρ 2 be smooth curves in P 2 (R m ) so that

ρ 1 (0) = ρ 2 (0) = µ.
(III.38)

P 2 (R m
) can be endowed with a Riemannian metric defined as follows [Otto 2001] :

∂ρ 1 ∂t (0), ∂ρ 2 ∂t (0) M,γ(t) = R m µ(x) v 1 (x), v 2 (x) dx, (III.39)
where v 1 and v 2 are two velocity fields in R m which satisfy the mass conservation equation III.32 w.r.t. ρ 1 and ρ 2 respectively, while guaranteeing that the associated instantaneous kinetic energies remain as small as possible [Villani 2003]. Then, the quadratic Wasserstein distance is precisely the associated Riemannian distance and most importantly, the unique geodesic between two distributions µ and ν is the previously defined curve χ * (t, .) (see Eq.III.30).

The process of interpolating two probability measures following this geodesic is known as the displacement interpolation because informally, the geodesic is obtained by advecting elementary mass particles from one distribution toward a target location defined by the optimal transport map in the other distribution, following straight lines [McCann 1997]. The simplicity of this mechanism makes it particularly appealing for numerical applications. Since W 2 is a geodesic distance, the curve χ * (t, .) solves the following optimization problem :

χ * (t, .) = argmin ω (1 -t)W 2 (µ, ω) 2 + tW 2 (ν, ω) 2 , (III.40)
∀t ∈ 0, 1 . Hence, χ * (t, .) can be formally interpreted as the barycenter of the probability measures µ and ν associated to the weights 1 -t and t in the Wasserstein metric. This notion has been extended to an arbitrary number of probability measures in [START_REF] Agueh | Martial Agueh et Guillaume Carlier. Barycenters in the Wasserstein space[END_REF] :

µ bar = argmin ω N i=1 w i W 2 (µ i , ω) 2 , (III.41)
with N i=1 w i = 1 and w i ≤ 0 ∀i ∈ 1, N . µ bar is referred as the (quadratic) Wasserstein barycenter of the probability measures (µ i ) 1≤i≤N and is precisely characterized in [START_REF] Agueh | Martial Agueh et Guillaume Carlier. Barycenters in the Wasserstein space[END_REF]. However, the choice of Wasserstein barycenters for averaging data samples can be motivated through the trivial example of Fig.VI.9. Indeed, the Wasserstein barycenter is somehow more meaningful than the multimodal euclidean barycenter. We postpone further discussions in this regard to Chapter VI. The weak gravitational lensing is one of the most promising tools to probe the dark matter distribution in the universe. The idea is to infer, from a billion images of galaxies, the shape distortions due to dark matter gravitational lensing and then estimate the dark matter mass density along different lines of sight. As mentioned in the introduction, the Euclid mission will provide the data for such a purpose. Nevertheless, galaxy images are distorted due to the PSF. Therefore, it is critical to know this distortion accurately. It can be modeled in first approximation as a convolution of the desired image by the PSF of the telescope, which is typically space and time-varying. In practice, isolated stars provide PSF measurements at different locations in the field of view. Nevertheless, these stars images can be aliased as it is the case in Euclid, given the CCD sensor sizes. On the other hand, the surveys are generally designed so that different images of the same stars are available and likely to be with different subpixel offsets on the sensor grid. Moreover, one may consider that nearby star images give to some extend the same local PSF. We, thus, can consider that different low-resolution versions of the same PSF are available in practice, so that one may apply an super-resolution (SR) method to recover aliased frequencies.

Chapitre IV

Multiple frame PSF super-resolution

This chapter precisely tackles this problem. The SR is a widely studied topic in general image processing literature. Yet, some methods have been specifically proposed for astronomical data. For instance, there is the software IMCOM [Rowe et al. 2011] and PSFEx, which proposes an SR option. The IMCOM provides an oversampled output image from multiple undersampled input images, assuming that the PSF is fully specified. Since it does not deal with the PSF restoration itself, we use PSFEx as our main reference. The PSFEx performs SR by minimizing the sum of two terms. The first term is a weighted quadratic distance, relating the underlying PSF to each of the low resolution measurements. The second term consists of the square l 2 norm of the difference between the underlying PSF and a smooth first guess. This term is meant for regularization. In the proposed algorithm, we introduce a new regularization scheme based on the optimization variable sparsity in a suitable dictionary.

Section A presents the general principle of SR along with some state-of-theart methods in astronomical domain. In Section B, we present the proposed algorithm in details, which is followed with some numerical experiments in Section C. We conclude by summarizing the main results and giving some perspectives.

A Super-resolution overview

A.1 Notations

The underlying high resolution (HR) image of size d 1 p 1 × d 2 p 2 is written in lexicographic order (for instance, lines after lines) as a vector of pixels values x = (x 1 , ..., x q ) T , where q = d 1 p 1 d 2 p 2 , and d 1 and d 2 are respectively the line and column downsampling factors. We consider n LR observations. The vector of pixels values y k = (y k1 , ..., y kp ) T denotes the k th LR observation written in lexicographical order with p = p 1 p 2 and k = 1...n.

A.2 Observation model

We assume that x does not change during the acquisition of the n LR images so that we have

y k = DB k M k x + n k , k = 1...n. (IV.1)
The variable M k is a warp matrix of size q × q. It represents the motions of the observations relative to each other and those in general need to be estimated. The variable B k is a blur matrix of size q×q. It accounts for different blurs (the atmosphere blur, which is particularly considerable for ground based telescopes, the system optics blur, the imaging system shaking etc.). The variable D is a matrix of size p × q, which simply realizes a downsampling operation. Finally, n k is a noise vector of q elements. This model is illustrated in In our case, we are interested in estimating the PSF, or in other terms, the telescope's contribution to the blur. We consider that the PSF varies slowly in the field so that the blocks "Warping" and "Blur" in Fig . IV.1 may be swapped, for slow motions between observations. Therefore, the block diagram can be adapted as in Fig. IV.2. This model still holds in presence of atmospheric blur (for ground-based telescopes) and jitter movements, if the LR images are extracted from the same exposure.

In the general case, the model (IV.1) may simply be written as

y k = W k x + n k , k = 1...n, (IV.2)
where W k is a p × q matrix accounting for warping, blur, and downsampling.

A.3 SR techniques in astronomy

Generally, SR techniques involve three steps, which may be combined or performed separately. The registration step consists in evaluating the relative motions between different observations, so that their samples can be arranged on a common grid. It is critical that the precision of this registration should be smaller than the pixels dimensions over the target upsampling factors. Since these relative motions are arbitrary, this grid is non-uniform. The next step would be to interpolate this grid in such a way to get a regularly sampled HR image. This image is blurry and noisy. Therefore, the final step is a restoration procedure. This scheme is summarized in Fig.

IV.3. In the next sections, we describe some SR techniques dedicated to astronomical images, but one may refer to [START_REF] Park | [END_REF]] for more details on various SR frameworks in general image processing literature.

A.3.1 Shift-and-add method

The most simple super-resolution method is certainly the shift-and-add method. It is performed in three steps. First, the images are upsampled to the target resolution. Then, they are shifted on a common grid and averaged. It has been used in astronomy for a long time, particularly for ground based telescopes. This method is simple, fast and is used for comparisons in the numerical experiments part. It has been shown that this method provides an optimal solution in the sense of the maximum likelihood (ML) with additive white Gaussian noise (WGN) and when only pure integer translation motions are considered with respect to the finer grid [Elad & hel Or 2001]. The interpolation operator should be D T (from Eq. IV.1), which comes down to a simple zero-padding interpolation. It has been shown in the same work that the matrix R = n k=1 M T k D T DM k is diagonal in this simple case. Thus, after registration and stacking of the interpolated images, each pixel value should be divided by the corresponding R diagonal coefficient.

A.3.2 PSFEx method

The software PSFEx is an open source program, which has been used in many projects such as the Dark Energy Survey [Mohr et al. 2012] or CFHTLS1 for PSF modeling. It takes a catalog of objects extracted from an astronomical image using SEXTRACTOR as input, which is also an open source tool for point sources (or stars) extraction. This catalog contains information about extracted point sources, such as SNR, luminosity, full width at half maximum (FWHM) , centroid coordinates, multiple flags related to saturation or blending, etc. Based on these measurements (performed in SEXTRACTOR) and some user provided parameters, PSFEx selects which sources are proper for PSF modeling. Afterwards, it constructs a PSF model, provides a fitting with an analytic function, and computes some of the PSF geometrical features. The PSF model construction may simply consist in optimally combining the input sources images for denoising, but it may also involve SR if these images are also undersampled. This SR functionality is our second reference for comparisons. These codes and the associated documentation may be found on the website http://www.astromatic.net/.

The desired HR image is now a matrix X of size dp × dp, where d is the downsampling factor and the k th observation Y k is a matrix of size p × p with k = 1...n. The coordinates of the centroid of the k th observation are denoted (i k , j k ). Assuming that the images are bandlimited, the samples of the LR images can be interpolated from the desired HR image, thanks to the Shannon sampling theorem [Shannon 1948]. In theory, this interpolation should involve a 2D sinus cardinal (sinc) kernel with an infinite support, which is not convenient for practical implementation. One can instead use a support compact function, which approximates the sinus cardinal. Let h(., .) denote such a function. The estimate of the sample (i, j) of the k th observation is given by

ŷk,ij = l m h [l -d(i -i k ), m -d(j -j k )] x ij . (IV.3)
Then, we can define the cost function

J 1 (X) = n k=1 p i=1 p j=1 (y k,ij -f k ŷk,ij ) 2 σ 2 k (IV.4)
where f k accounts for possible luminosity differences. The parameter σ 2 k is related to the local background variance and any other uncertainty on the pixel value. These parameters need to be estimated.

The PSFEx uses a Lanczos4 interpolant, which is defined in 1D as

h 1D (x) =    1 if x = 0 sinc(x) sinc(x/4) if 0 < |x| < 4 0 else, (IV.5) so that h(x, y) = h 1D (x)h 1D (y).
Finally, PSFEx minimizes the cost function defined as

J 2 (∆) = J 1 (∆ + X (0) ) + λ ∆ 2 2 , (IV.6)
where ∆ = X-X (0) and X (0) is a median image computed from the LR observations. This second term is meant for regularization purposes if the problem is ill-conditioned or undetermined (n < d 2 ). One particularity of point-source images (see Fig.

IV.4) is that one knows that the light comes from a single point, which is generally inferred to be the light blob's centroid. Thus, one only needs the images centroid coordinates to perform the registration, which is implicitly done in Eq. IV.3. where x and ŷk are, respectively, the desired matrix and the k th observation estimate written this time as column vectors in lexicographic order (we used the lines order) ; besides, H k is a Toeplitz matrix [START_REF] Gray | [END_REF]] of size d 2 p 2 ×d 2 p 2 , which contains the values of the kernel h(., .) appearing in Eq. IV.3 and D is a decimation matrix of size p 2 × d 2 p 2 with d being the downsampling factor. We are assumming that the images are squared for convenience but they might be non-squared as well. In the same way, we can redefine the objective function of Equation IV.4 as

J 1 (x) = 1 2 n k=1 y k -f k DH k x 2 2 /σ 2 k , (IV.8)
where y k is the k th observation rewritten consistently with x and ŷk . The function J 1 (x) is nothing but the negative log-likelihood associated with the observation model in the case of an uncorrelated Gaussian noise that is stationary for each observation up to a positive scalar factor. It can be written in an even more compact way as

J 1 (x) = 1 2 Σ -1 y -Σ -1 FWx 2 2 , (IV.9)
where W is obtained by concatening vertically the matrices DH k , F is a diagonal matrix constructed by repeating the coefficients f k p 2 times for k = 1...n, and Σ is constructed the same way using the coefficients σ k . Therefore, we can simply write

J 1 (x) = 1 2 z -Mx 2 2 , (IV.10)
where M is a matrix of size np 2 × d 2 p 2 . The SPRITE method constrains the minimization of this objective function using an analysis prior. Instead of using a single lagrangian multiplier as in Eq. II.41, the analysis coefficients has individual weights κλ i , where i is the coefficient index in the transform domain. This leads to the following formulation of the problem :

min ∆ J 1 (∆ + x (0) ) + κ λ Φ∆ 1 , (IV.11)
where ∆ = xx (0) is defined as in IV.6, x (0) is a first guess, and λ is now a vector of the same size as Φ∆, denoting the pointwise product.

Additionally, the PSF or equivalently the telescope optical impulse response is by definition a positive valued function. Therefore, we want the reconstruction x to have positive entries. This additional constraint is integrated as follows :

min ∆ J 1 (∆ + x (0) ) + κ λ Φ∆ 1 s.t.∆ ≥ -x (0) , (IV.12)
where ≥ is a pointwise inequality, or equivalently,

min ∆ J 1 (∆ + x (0) ) + κ λ Φ∆ 1 + ι S x (0) (∆), (IV.13)
where S x (0) is the set of vectors t ∈ R d 2 p 2 satisfying the pointwise inequality t ≥ -x (0) and ι S x (0) is its indicator function (see Appendix B). The impact of this constraint is emphasized in Appendix C. As we show in the Section B.3.3, the choice of the parameters λ and κ relies on the noise expected on the analysis coefficients of the solution estimate. The choice of a vector regularization parameter rather than a single scalar is precisely motivated by the fact that this noise might be non-stationary.

As stated before, the problem IV.11 is not equivalent to its synthesis version, if the dictionary Φ is redundant. The synthesis prior is expected to be efficient if the desired solution can be accurately written as a sparse linear combination of the chosen dictionary atoms, which we cannot assume to be true for every PSF profiles. In contrast, the analysis prior appears to be more flexible. Moreover, in the cases where the problem would be ill-conditioned or underdetermined, the analysis prior would definitely be more suitable since it involves far less variables.

A similar 1 penalty has already been applied for SR. An example may be found in [Yamagishi et al. 2012], where the cost function is minimized using variants of the alternating direction method of multipliers (ADMM). Moreover, advantages of such approaches over quadratic regularizers have been shown in many related problems in image and signal processing.

The use of the 1 norm as a relaxation for an l 0 penalty has a well-known drawback, which is that it tends to bias the solution. Indeed with a 1 norm penalty, the problem resolution involves soft thresholding operations, which affect both weak and strong entries unlike a hard thresholding, which would only affect the weak and therefore unwanted entries. Formal definitions of soft and hard thresholding are given in Appendix B. This is particularly unsuitable for scientific data analysis. The reweighting 1 minimization proposed in [Candès et al. 2008] is one way to tackle this issue, while staying in the proof of convergence sets. Indeed, it consists of solving a succession of 1 minimization problems of the form

min ∆ J 1 (∆ + x (0) ) + κ w (k) λ Φ∆ 1 + ι S x (0) (∆), (IV.14)
where w (k) is a weighting vector for the transform coefficients at the k th minimization. Each entry of w (k) is calculated as a decreasing function of the corresponding transform coefficient magnitude in the (k -1) th minimization. This way, the strong transform coefficients are less penalized than the weaker ones in the new minimization. One may refer to Appendix C for quantitative study of the reweighting effect.

B.2 Algorithm

In the proposed method, the reweighting scheme is performed according to [Candès et al. 2008] :

1. Set k = 0, for each entry in the weighting vector w (k) , set w

(k) j = 1.
2. Solve the problem IV.14 yielding a solution ∆ (k) .

Compute α

(k) = Φ∆ (k) .
4. Update the weight vector according to w

(k+1) j = 1 1+|α (k) j |/3σ j
, where σ j is the noise standard deviation expected at the j th transform coefficient, see section B.3.3. 5. Terminated on convergence or when reaching the maximum number of iterations ; otherwise, go to step 2.

The step 2 resolution is detailed in Algorithm IV.1. We use the generalized forward-backward splitting introduced in [Raguet et al. 2011]. It requires the computation of proximity operators associated with the regularization functions in IV.14. One may refer to Appendix B for an introduction to proximal calculus.

Since the dictionary Φ is redundant, we do not have a closed-form expression for the proximity operator. Yet, it can be calculated as

prox µ ω 1 κ Diag (w (k) λ)Φ 1 (x) = x -Φ T û û = argmin |u j |< µ ω 1 κw (k) j λ j 1 2 x -Φ T u 2 2 (IV.15)
where û can be estimated using a forward-backward algorithm [START_REF] Bauschke | [END_REF] as follows :

1. Set p = 0, initialize u 0 = 0. 2. ũp+1 = u p + µ prox Φ x -Φ T u p . 3. u p+1 = ũp+1 -SoftThresh µ ω 1 κw (k) λ ũp+1 .
4. Terminate on convergence or when reaching the maximum number of iterations, otherwise go to step 2.

Algorithme IV.1 Weighted analysis-based ∆ (k) recovery

Paramètre(s) d'entrée :

A first guess estimate of the super-resolved image x (0) . A weight vector w (k) . Sparsity constraint parameter κ.

A dictionary Φ.

Auxiliary variables z 10 , z 20 ∈ R d 2 p 2 . ω 1 , ω 2 ∈]0, 1[ s.t. ω 1 + ω 2 = 1, λ > 0 (see Section B.3.4). 1: Initialize d 0 = ω 1 z 10 + ω 2 z 20 . 2: Pour n = 0 to N max -1 faire 3: z 1n+1 = z 1n + λ(prox µ ω 1 κ Diag (w k λ)Φ 1 (2d n -z 1n -µ∇J 1 (d n + x (0) ) - d n ) 4: z 2n+1 = z 2n + λ(prox µ ω 2 ι S x (0) (∆) (2d n -z 2n -µ∇J 1 (d n + x (0) ) -d n ) 5: d n+1 = ω 1 z 1n+1 + ω 2 z 2n+1 6: Return : ∆ (k) = d Nmax .
The thresholding operator SoftThresh is defined in Appendix B. The operator prox

µ ω 2 ι S x (0)
is simply the orthogonal projector onto the set S x (0) defined in the previous section ; it is given explicitly Appendix B.

A full description of SPRITE is provided in Algorithm IV.2.

Algorithme

∆ (k) = argmin ∆ J 1 (∆ + x (0) ) + κ w (k) λ Φ∆ 1 + ι S x (0) (∆) (see Algorithm IV.1) 7: α (k) = Φ∆ (k) 8: w (k+1) j = 1 1+|α (k) j |/3σ j 9: Return : x = ∆ (Kmax-1) + x (0) .

B.3 Parameter estimation

The data fidelity term in the problem IV.14 is defined as

J 1 (x) = 1 2 n k=1 p i=1 p j=1 (y k,ij -f k ŷk,ij ) 2 σ 2 k ŷk,ij = l m h [l -d(i -i k ), m -d(j -j k )] x ij , (IV.16)
where x = (x ij ) 1≤i,j≤dp is the desired image and h(., .) is a 2D Lanczos kernel.

As stated in Section B.1, it can be written as

J 1 (x) = 1 2 z -Mx 2 2 , (IV.17)
if we write x in lexicographic order as a vector. The following parameters are required :

-for the data fidelity term parameters the noise standard deviation in the LR images σ k , the photometric flux f k , and the shift parameters (i k , j k ) ; -a first guess x (0) (see problem IV.14) ;

-the sparsity constraint parameters κ, λ = (λ i ) i and the dictionary Φ ;

-algorithmic parameter such as the gradient step size µ, the relaxation parameter λ in Algorithm IV.1, and the gradient step size µ prox in Equation IV.15 resolution.

B.3.1 Data fidelity parameters

Noise standard deviations At the first step of the algorithm, the noise standard deviations in the low resolution images can be robustly estimated using the median absolute deviation (MAD) estimator [START_REF] Starck | [END_REF].

Subpixel shifts

The subpixel shifts between the images are estimated based on the low resolution images centroids positions. Those are calculated on the low resolution images after a hard thresholding operation. For the image x i , the threshold is chosen as

k = min(4σ i , max(|x i |) σ i -1 σ i ). (IV.18)
In this way, we only keep pixels with a high SNR for the centroid estimation. We then estimate the centroid positions using the iteratively weighted algorithm introduced in [START_REF] Baker | [END_REF]. The thresholding operation undoubtedly biases the estimated centroid position, but the resulting estimated shifts are expected to be unbiased, up to the finite sampling and noise effects.

Photometric flux

The flux parameters are calculated by integrating the low resolution images on a fixed circular aperture centered on their centroids estimates. At Euclid resolution (see Section C), we obtained quite accurate flux estimates in simulations using a radius of 3 pixels for the aperture. These parameters define the matrix M in Eq. IV.17.

All these parameters are automatically calculated without requiring any user input.

B.3.2 First guess computation

As one can see in the Algorithm IV.2, the final image is computed as

x = ∆ (Kmax-1) + x (0) . (IV.19)
This implies that the noise and any artifact in x (0) which does not have a sparse decomposition in Φ will be present in the final solution. Therefore, one has to be careful at this step. To do so, we compute a noisy first guess x (0) n using a shift-and-add, as presented in Section A.3.1. Then we apply a wavelet denoising to x

n . In other terms, we transform x (0) n in a "sparsifying" wavelet dictionary W. We threshold each wavelet scale in such a way to keep only the coefficients above the noise level expected in the scale. Finally, we apply a reconstruction operator, which is a dictionary W verifying WW = I q to the thresholded coefficients (see [START_REF] Starck | [END_REF]). We note β = (β i ) i , a vector made of the denoising thresholds for each wavelet scale. We set β i at 5σ i , where σ i is the noise standard deviation in the i th wavelet scale. The first guess is finally computed as

x (0) = W HardThresh β Wx (0)
n , (IV.20) which robustly removes the noise without breaking important features. One can refer to Appendix C for (σ i ) i estimation.

B.3.3 The choice of dictionary and regularization parameter

Regularization parameter The regularization parameter κ can be set, according to a desired level of significance. Indeed, it can be seen that the transform domain vector û is constrained into weighted l ∞ -ball of radius µκ in Equation IV.15 and can be interpreted as the non-significant part of the wanted signal current estimate. To set this radius according to the expected level of noise for each transform coefficient, we propagate the noise on the data vector z from Equation IV.17 through µΦM T M and estimate its standard deviation at each transform coefficient, which sets the parameters λ j . In practice, this can be done in two ways. We can either run a Monte-Carlo simulation of the noise in z and take the empirical variance of the sets of realizations of each transform coefficient. On the other hand, if Φ is a wavelet dictionary and if the noise is expected to be stationary in each wavelet scale, then we only need to compute a single standard deviation per scale. This can be done by estimating the noise in each scale of the wavelet transform of the gradient at each iteration (up to the factor µ) using a MAD, for instance. Indeed, the residual z -M(d n + x (0) ) tend to be consistent with the noise in z, so that it can be used as a noise realization. With a stationary noise in each input image, the two approaches give very close estimates of the noise standard deviation and the second one is far less demanding in terms of complexity. As a result, coefficients below κλ j are considered as part of the noise and one only needs to set the global parameter κ to tune the sparsity constraint according to the noise level.

Dictionary The choice of the dictionary impacts the performance of the algorithm. We considered two transforms : a biorthogonal undecimated wavelet transform with a 7/9 filter bank and the second generation starlet transform [Starck et al. 2011]. These two transforms are generic and not specifically tuned to a given PSF profile.

B.3.4 Algorithmic parameters

Gradient steps sizes The gradient step size µ in Algorithm IV.1 needs to be chosen just in 0, 2/ρ(M T M) , where ρ(.) denotes the spectral radius of a square matrix. In the same way, µ prox needs to be chosen in 0, 2/ρ(ΦΦ T ) .

Relaxation parameter

The parameter λ in Algorithm IV.1 needs to be chosen in 0, min 3 2 , 1+2/ρ(M T M)µ 2 [Raguet et al. 2011]. This parameter tunes the updating speed of the auxiliary variables in Algorithm IV.1. In practice, we use µ = 1/ρ(M T M) and λ = 1.4.

B.3.5 User parameters

It is important to mention that the user only has to set the parameter κ and the dictionary with the other parameters being automatically estimated. In all our experiments, we took κ = 4, which is quite convenient, if we assume Gaussian noise in the data. The dictionary choice will be emphasized in the next section.

C Numerical experiments

This section presents the data used, the numerical experiments realized as a mean to compare three SR techniques (shift-and-add, PSFEx method, and our method) and the results.

C.1 Data set

The PSFs provided are optical PSFs computed using a fast Fourier transform of the exit pupil. They are not a system PSF, so they do not include a jitter or detector response. A set of PSFs covering the whole field of view is provided. They are monochromatic PSFs at 800nm and are derived from tolerance analysis. They account for manufacturing and alignments errors and thermal stability of the telescope. Manufacturing and alignment errors are partially compensated by a best focus optimization, while thermal stability effects are simulated by a small displacement of the optics that are not compensated on a short-time scale. The optical model used is dated from 2011 and is prior to the current reference model (provided by Astrium, which has been awarded the payload module contract in 2013). In particular, the 2011 model does not contain the latest definition of the pupil mask. The pupil, however, includes central obscuration and a three-vane spider. This is the model that has been used for the science feasibility studies that led to the acceptance of the Euclid mission.

C.2 Simulation

In the Euclid mission, the actual sampling frequency is about 0.688 times the Nyquist frequency that we define as twice the telescope spatial cut-off frequency [Cropper 2013]. Therefore we target an upsampling factor of 2, which gives a sufficient bandpass to recover the high frequencies. The PSF is typically space-varying, and this is particularly true for wide field of view instruments as Euclid telescope [ESA/SRE 2011]. Thus, the data set contains simulated PSF measurements on a regular 18 × 18 grid on the field of view. The original PSF models are downsampled to twice Euclid resolution. For each PSF, four randomly shifted "copies" are generated and downsampled to Euclid resolution (see Fig. We define the signal level as its empirical variance that is calculated in a 50 × 50 patch centered on the HR PSF main lobe. For each algorithm, we used these four images to reconstruct a PSF which has twice their resolutions in lines and columns.

C.3 Quality criterion

For an image X = (x ij ) i,j , the weighted central moments are defined as

µ p,q (X) = i j (i -i c ) p (j -j c ) p f ij x ij (IV.21)
with (p, q) ∈ N2 , (i c , j c ) are the weighted image centroid coordinates, and F = (f ij ) i,j is an appropriate weighting function (typically a Gaussian function). The ellipticity parameters are then defined as follows :

e 1 (X) = µ 2,0 (X) -µ 0,2 (X) µ 2,0 (X) + µ 0,2 (X) (IV.22) e 2 (X) = 2µ 1,1 (X) µ 2,0 (X) + µ 0,2 (X)
.

(IV.23)

The vector ε = [e 1 , e 2 ] is an important tool, since if measured on a large set of galaxies, it can be statistically related to the dark matter induced geometrical distortions and finally its mass density. Furthermore, this ellipticity parameters are magnitude invariant and approximately shift invariant. The error on ellipticity is therefore an interesting criteria for quality assessment. Thus, we used the mean absolute error for each ellipticity parameter,

E j = 1 n n i=1 |e j (X i ) -e j ( Xi )|, j = 1, 2 (IV.24)
and the associated empirical standard deviations. Moreover, the PSF size is also an important characteristic of the PSF kernel. For example, it has been shown in [START_REF] Paulin-Henriksson | [END_REF]] that the PSF size largely contributes to the systematic error in weak gravitational lensing surveys. Therefore, we use it in quality assessment by computing the mean absolute error on the full width at half maximum (FWHM). The FWHM is estimated by fitting a modified Lorentzian function on the PSF images. We used routines from a publicly available library 2 .

C.4 Results and discussion

The Figures IV.5 and IV.6 show a simulated PSF that is sampled at almost the Nyquist rate and the two LR shifted and noisy PSF derive, with SNR of around 30dB. The Figure IV.7 shows an example of super-resolved PSF at 30dB of SNR, from four LR images and the corresponding error maps that are defined as the absolute value of the difference between the original high resolution noise free PSF and the PSF reconstructions for each algorithm. This error map standard deviation is at least 30% lower with SPRITE.

Figure IV.8 shows smaller errors and errors dispersions are achieved with SPRITE algorithm, especially at low SNR. One can note that the dispersion is slightly smaller with biorthogonal undecimated wavelet. The error on the FWHM given in percent on Fig. IV.9 is smaller with SPRITE. In practice, there is more variability in the PSF (wavelength and spatial dependency, time variations...) so that the real problem will be more underdetermined. Thanks to the multiple exposures on the one hand, and that the spatial variations of the PSF are expected to be slow on the other hand, the real problem could actually be very well constrained. Moreover, these results suggest that even better results could be achieved by using more adapted dictionaries, built either from PSF model or through a dictionary learning algorithm [Beckouche et al. 2013a].

D Complexity and performances

The simulations were run on a typical desktop computer. Let us suppose that we have n LR images of sizes p 1 × p 2 and that we choose an upsampling factor d in lines and columns. As stated before, we took p 1 = p 2 = 84, n = 4, and d = 2 in our numerical experiments. Under this setting, it takes roughly 60s and 1GB of physical memory to compute a super-resolved PSF. More generally, the computational complexity of the algorithm is in O(np 1 p 2 d 2 log(p 1 p 2 d 2 )), which is related to the implementation of the matrices M from Eq. IV.10 and M T using FFT.

E Software

Following the philosophy of reproducible research [START_REF] Buckheit | WaveLab and Reproducible Research[END_REF], the algorithm introduced in this paper and the data used are available at http://www.cosmostat.org/sprite.html. We used the following calls for the SPRITE executable :

-run_sprite -t 2 -s 4 -r 2 -F -N data_file output_file out-put_directory for the second genration Starlet transform ;

-run_sprite -t 24 -s 4 -r 2 -F -N data_file output_file out-put_directory for the undecimated biorthogonal wavelet transform. The options "-s" and "-r" set the parameter κ (see Section B.3.3) and the upsampling factor for both lines and columns respectively. The options "-F" and "-N" indicate that the photometric flux and the noise might have different levels in the LR images and need to be estimated.

F Conclusion

We introduced SPRITE, which is a super-resolution algorithm based on sparse regularization. We show that adding a sparse penalty in the recovery leads to far better accuracy in terms of ellipticity error, especially at low SNR. Quantitatively, we achieved -a 30% lower error on the reconstruction itself at 30dB of SNR ;

-around 6dB less than other methods on the shape parameters, which corresponds to a factor of e 6 on a linear scale, at 10dB of SNR ;

-6% of error on the FWHM in average, which 2% less than PSFEX.

However this algorithm does not handle the PSF spatial variations. Thus one natural extension of this work would be to simultaneously perform superresolution and dimensionality reduction assuming only one LR version of each PSF, as it is the case in practice strictly speaking, but by using a large PSF set. This is precisely the subject of the next chapter. In the previous chapter, the PSFs super-resolution (SR) where performed, assuming that multiple measurements of the same PSF were available. However in practice, one does not have such multiple measurements strictly speaking, because the PSF is space and time variant, as emphasized in Chapter I.

In this chapter, we consider precisely a setting where the PSF is space variant and we want to get an accurate modeling at high resolution of the PSF field, assuming we have under-sampled measurements of different PSFs in the observed field. We assume that the PSFs vary slowly across the field. Intuitively, this implies a compressibility of the PSFs field, which leads us to the question of what would be a concise and easily understandable representation of a spatially indexed set of PSFs.

A The PSF field model

A.1 The observation model

We assume that we have an image I, which contains p unresolved objects such as stars, which can be used to estimate the PSFs field. Noting y k one of these p objects at spatial position u k , y k is therefore a small patch of I with n y pixels, around the spatial position u k . We will write y k as a 1D vector. The relation between the "true" PSF x k and the noisy y k observation is

y k = M k x k + n k (V.1)
where M k is a linear operator and n k is a noise that we assume to be Gaussian and white. We will consider two kinds of operators in this paper : the first one is the simple case where M k = I nx and we have the number of pixels n x in x k is equal to n y , and the second one is a shift+downsampling degradation operator and n x = m 2 d n y , where m d is the downsampling factor in lines and columns, with m d ≥ 1.

Noting Y = [y 1 . . . y p ] the matrix of n y lines and p columns of all observed patches, X = [x 1 . . . x p ] the matrix n x ×p of all unknown PSFs, we can rewrite Eq. V.1 as

Y = F(X) + N (V.2)
where

F(X) = [M 1 x 1 , . . . , M p x p ].
This rewriting is useful because, as we discuss in the following, the different PSFs x k are not independent, which means that the problems of Eq. V.1 should not be solved independently for each k. In other terms, the vectors (x k ) 1≤k≤p belong to a specific unknown manifold that needs to be learned by using the data globally.

A.2 The data model

Let Ω be a r dimensional subspace of R nx embedding the PSFs field. We assume that there exists a continuous function f : E → Ω, so that f (u k ) = x k , ∀k ∈ 1, p . The regularity of f translates the correlation of the data in space (and time).

Let (s i ) 1≤i≤r be a basis of Ω. By definition, we can write each x k as a linear combination of the s i , x k = r i=1 a ik s i , k = 1 . . . p, or equivalently

X = SA (V.3)
where S = [s 1 , . . . , s r ] and A is a r×p matrix containing the coefficients A[:, k] of the vectors x k (k = 1 . . . p) in the dictionary S. Each column of the matrix S, that we also refer to as an atom, can be seen as an eigen PSF, i.e. a given PSF's feature distributed across the field.

A.3 The inverse problem

We need therefore to minimize Y-F(X) 2 F , which is an ill posed problem due to both the noise and the operator F, . F denoting the Frobenius norm of a matrix. There are several constraints that may be interesting to use in order to properly regularize this inverse problem :

-positivity constraint : the PSF x k should be positive ; -low rank constraint : as described above, we can assume that x k = r i=1 a ik s i , which means that we can instead minimize

min A,S Y -F(SA) 2 F ; (V.4)
we assume that r min(n, p) ; this dimension reduction has the advantage that there are much less unknown to find, leading to more robustness, but the problem is now that the cost function is not convex anymore ;

-smoothness constraint : we can assume that the vectors x k are structured ; the low rank constraint does not necessarily impose x k to be smooth or piece-wise smooth ; adding an additional constraint on S atoms, such as a sparsity constraint, allows to capture spatial correlations within the PSFs themselves ; an additional dictionary Φ s can therefore be introduced which is assumed to give a sparse representation of the vectors s k ;

-proximity constraint : we can assume that a given x k at a position u k is very close to another PSF x k at position u k if the distance between u k and u k is small ; this means that the field f must be regular ; this regularity can be forced by adding constraints on the lines of the matrix A ; indeed, the p values relative to a line A[i, :] correspond to the contribution of the ith eigen PSF to locations relative to the spatial positions U.

We show in section C how these four constraints can be jointly used to derive the solution. Let first review existing methods susceptible to solve this problem.

B Related work

In all this section, Y refers to the observed data set (y k ) 1≤k≤p . In the first part, the aforementioned degradation operator F is simply the identity. Therefore we review some dimension reduction methods. In the second part F is a shifting and downsampling operator ; we present a PSF modeling software dealing with this more constraining setting.

B.1 Dimension reduction

The principal components analysis is certainly one of the most popular mathematical procedure in multivariate data analysis and especially, dimension reduction. In our case, we want to represent Y's elements using r vectors, with r ≤ max(p, n y ). A PCA gives an orthonormal family of r vectors in R ny so that the total variance of Y along these vectors directions is maximized. By definition, the PCA looks for redundant features over the whole data set. Therefore, in general, the principal components neither capture localized features (in sense of E) nor have a simple physical interpretation.

In [START_REF] Wang | [END_REF], a "regularized" PCA is proposed to address this shortcoming for spatial data analysis in atmospheric and earth science. Indeed, as a PCA, the method solves the following problem,

min A Y -YA T A 2 F , s. t. AA T = I r , (V.5)
for some chosen small r. Moreover, it jointly imposes a sparsity constraint and a smoothing penalties with respect to the space E, on the matrix A lines. This way, with the right balance between those two penalties, one favors the extraction of localized spatial features, making the interpretation of the optimal A easy. Yet, there is no obvious way of setting the sparsity and smoothness parameters, which are crucial ; moreover, unless the data actually contain spatially localized and non-overlapping features, the coupled orthogonality and sparsity constraint is likely to yield a biased approximation of the data.

In the context of remote sensing and multi-channel imaging, two ways of integrating spatial information into PCA are proposed in [Cheng 2006] ; the set Y is made of multi-channel pixels. In the first way, the author introduces a weighting matrix indicating the relative importance of each pixel. For instance, the weight of a given pixel can be related to its distance to some location of interest in E. Then, the computation of the covariance matrix of image bands is slightly modified to integrate this weighting. This idea is close to the methodology proposed in [Harris et al. 2011]. As a consequence, one expects to recover spectral features spatially related to some location of interest within the most important "eigen-pixels". Yet, we do not have any specific location of interest in E and we rather want to recover relevant features across the whole data set.

The second approach aims at taking into account the spatial associations and structural properties of the image. To do so, modified versions of the image bands covariance matrices are calculated, with increasing shifts between the bands, up to a predetermined maximum shifting amplitude. These covariance matrices, including the "regular" one, are averaged and the principal components are finally derived. Intuitively, one expects the spectral features present in structured images regions to be strengthened and therefore captured into the principal components. However, we consider a general setting where the data are randomly distributed with respect to E, which makes the shifted covariances matrices ill-defined.

A review of PCA applications and modifications for spatial data analysis can be found in [START_REF] Demšar | [END_REF].

In case the data lie on or are close to a manifold M of dimension r embedded in R n , one can consider using one of the numerous non-linear dimension reduction algorithms published in the manifold learning literature, such as GMRA [Allard et al. 2012], [START_REF] Maggioni | [END_REF]]. The idea is to partition the data in smaller subsets of sample close to each other in the sense of the manifold geometry. From this partionning, the manifold tangent spaces are estimated at subsets locations ; estimates are then simply given by the best regressions of these subsets with r-dimensional affine subspaces. The method includes some multiresolution considerations that are not relevant to our problem. This procedure provides a dictionary in which each of the original samples need at most r elements to be represented. Moreover, the local processing of the data, which is necessary in this setting because of the manifold curvature, makes this approach somehow compatible with the considered problem. Indeed, by hypothesis, the closer two samples will be in sense of E, the closer they will be in R n , and the more likely they will fall into the same local cluster.

Another interesting alternative to the PCA can be found in [START_REF] Lee | [END_REF]. This construction called "Treelets" extracts features by uncovering correlated subsets of variables across the data samples. It is particularly useful when the sample size is by far smaller than the data dimensionality (p n y ), which does not hold in the application we consider in the following.

B.2 Super-resolution

In this subsection, F takes the following form :

F(X) = [M 1 x (c) 1 , . . . , M p x (c) p ], (V.6)
where M i is a warping and downsampling matrix. Since we consider a set of compact objects images, the only geometric transformation one has to deal with for registration is the images shifts with respect to the finest pixel grid, which can be estimated using the images centroids as shown in Chapter IV.

To the best of our knowledge, the only method dealing with this specific setting is the one used in the PSF modeling software PSFEx [Bertin 2011]. This method solves a problem of the form :

min ∆ S 1 2 Y -F((∆ S + S 0 )A) 2 F + λ ∆ S 2 F . (V.7)
S 0 is a rough first guess of the model components. Each line of the weight matrix A is assumed to follow a monomial law of some given field's parameters. The number of components is determined by the maximal degree of the monomials. For instance, let say that we want to model the PSFs variations as a function of their position in the field with monomials with degrees up to 3, then :

-one needs 6 components corresponding to the monomials 1, X, X 2 , Y, XY and Y 2 ;

-assuming that the i th PSF in Y's columns order is located at u i = (u ix , u iy ) then the i th column of A is given by a

(c) i = [1, u ix , u 2 ix , u iy , u ix u iy , u 2 iy ]
T up to a scaling factor. This method is used for comparisons in the Numerical experiments part.

C Resolved Components Analysis

C.1 Matrix factorization

We have seen that we can describe the PSFs field f as

[f (u 1 ), . . . , f (u p )] = X = SA. (V.8)
The matrix S is independent of the spatial location, and the i th line of A gives the contribution of the vector s i to each of the samples. As discussed in section A.3, the field's regularity can be taken into account by introducing a structuring of the matrix A. We can write : Figure V.1 -Data matrix factorization : the j th sample, which is stored in the j th column of X is linear combination of S columns using A's j th column coefficients as the weights ; similarly, the j th line of A is a linear combination V T 's lines, using α's j th line coefficients as the weights.

A[i, :] T = N l=1 α il υ l , i = 1 . . .

Physical interpretation

An interesting way to well interpret A is to consider the ideal case where the measurements are distributed following a regular grid of locations U. In this case, we can expand the vector A[i, :] T using the Discrete Cosine Transform (DCT), and vectors υ i in Eq. V.9 are regular cosine atoms, and the column index of the matrix is related the frequency. Hence, lines relative to high frequencies will be related to quicky varying PSF components in the field, while lines related to low frequencies will be related to PSFs stable components. In practice, the sampling is not regular and the DCT cannot be used, and V has to be learned in a way to keep the harmonic interpretation valid. We want some lines A[i, :] to describe stable PSFs components on the FOV, and other to be more related to local behavior.

C.2 The proximity constraint on A

As previously mentioned, we want to account for the PSFs field's regularity by constraining A's lines. Specifically, we want some lines to determine the distribution of stable features across the PSFs field while we want other lines to be related to more localized features. In order to build this constraint, let first consider the simple case of a one dimensional field of regularly spaced PSFs.

C.2.1 Regularly distributed observations

We first assume that E = R. We suppose that p = 2k + 1, for some integer k and we consider the 1D vector ψ e,a = (ψ i ) 1≤i≤p defined as follows :

ψ i = ψ p-i+1 = -1/|u i -u k+1 | e if i = k + 1, (V.10) ψ i = p j=1 j =k+1 a/|u j -u k+1 | e , otherwise (V.11)
for some positive reals e and a. We suppose that ψ e,a is normalized in l 2 norm. We refer to this family of signals, parametrized by e and a as "notch filters", in reason of their frequency responses shapes. Some examples can be found in Fig.

V.2b. One can observe that ψ 1,1 is essentially a high pass filter. As e increases, the notch structure clearly appears, with an increasing notch frequency. It is clear that, for a vector v, minimizing the functional

Ψ e,a (v) = v ψ e,a 2
2 promotes vectors with spectra concentrated around the notch frequency corresponding to the chosen values of e and a. We can directly use this family of filters to constraint A as follows : we define the functional

Ψ : M rp (R) → R + , A → r i=1 Ψ e i ,a (A[i, :]), (V.12)
where (e i ) i is a set of reals verifying 0 ≤ e 1 < e 2 < • • • < e r and a ∈ [0, 2[. Because the notch frequency increases with e i , minimizing Ψ promotes varying level of smoothness of A's lines, which is what we wanted to achieve. The filter ψ e,a and the functional definitions can be extended to higher dimensions for some positive λ. Yet, in practical applications, the observations are in gene-ral irregularly distributed. In the next section, we propose a slightly different penalty which is usable for arbitrary observations distributions.

C.2.2 General setting

Let define the functional

Ψ e,a : R p → R + , v → p k=1 ( p i=1 i =k av k -v i u k -u i e 2 ) 2 , (V.14)
where e and a are positive reals. Minimizing Ψ e,a (v) tends to enforce the similarity of close features, with respect to E ; in other terms, the more u k -u i 2 is large, the less important is (av k -v i ) in Ψ e,a (v) and e somehow determines the radius of similarity. For e > 1, Ψ e,a ≈ Ψ e,a because of the uniform spacing of the values u i and the decay of

1 u k -u i e 2
, for sufficiently high p ; we give more details on this approximation in E in the 1D case. However unlike Ψ e,a , the functional Ψ e,a is still relevant without the uniform sampling hypothesis and we expect qualitatively the same behavior as Ψ e,a with respect to the frequency domain if the data sampling is sufficiently dense. Therefore we use Ψ e,a instead of Ψ e,a in the functional Ψ of Eq.V.13. Besides, we use the term frequency even for randomly distributed samples.

C.2.3 Flexible penalization : the redundant frequencies dictionary V

The efficiency of the regularization of the problem V.13 relies on a good choice of the parameters e 1 , . . . , e r and a. Indeed, if the associated notch frequencies does not match with the data set frequency content, the regularization will more or less bias the PSF estimation depending on the Lagrange multiplier λ. Besides, setting this parameter might be tricky. We propose an alternate strategy for constraining A, which leads to the factorization model introduced in Section C.1 and still builds over the idea of notch filters.

For v ∈ R p we can write

Ψ e,a (v) = P e,a v 2 2 , (V.15)
where P e,a is a p × p matrix defined by

P e,a [i, j] = - 1 u i -u j e 2 if i = j, (V.16) P e,a [i, i] = p j=1 j =i a u i -u j e 2 ,
(V.17)

(i, j) ∈ 1, p 2 . Therefore, Ψ e,a (v) = v T Q e,a v, (V.18)
where Q e,a = P T e,a P e,a and is symmetric and positive. We consider the singular values decomposition (SVD) of Q e,a : Q e,a = V e,a D e,a V T e,a . The diagonal values of D e,a are sorted in decreasing order. We note d e,a the vector made of these diagonal values, so that

d e,a [1] ≥ • • • ≥ d e,a [p] ≥ 0. Consi- dering the reduced form Ψ e,a (v) = p i=1 d e,a [i] v, V e,a [:, i] 2
, it is clear that minimizing Ψ e,a (v) promotes vectors correlated with Q e,a last eigenvectors. In the case of regular sampling with respect to E, these eigenvectors are the harmonics close to the notch frequency of ψ e,a . We can rewrite the functional Ψ accordingly :

Ψ(A) = r i=1 p j=1 d e i ,a [j] v, V e i ,a [:, j] 2 . (V.19)
It is clear from this expression that minimizing Ψ(A) enforces the selection of the eigenvectors associated with the lowest eigenvalues in the set (d e i ,a [j]) i,j for describing A's lines. This can be seen as a sparsity constraint over A's lines with respect to the atoms (V e i ,a [:, j]) i,j ; yet, the small subset of atoms which will carry most of the information is somehow predefined through the eigenvalues (d e i ,a [j]) i,j . This is unsuitable if the notch filters parameters are poorly selected ; on the contrary, one would like to select in a flexible way the atoms which fit the best the data.

Let suppose that we have determined a set of parameters (e i , a i ) 1≤i≤r so that the filters ψ e i ,a i notch frequencies would cover the range of significant frequencies (with respect to E) present in the data. As previously, we note (V e i ,a i ) 1≤i≤r the eigenvector's matrices associated with the operators Ψ e i ,a i . We note V = [V e 1 ,a 1 , . . . , V er,ar ]. Considering the preceding remark, we introduce the following problem :

min α,S 1 2 Y -F(SαV T ) 2 F s.t. α[l, :] 0 ≤ η l , l = 1 . . . r (V.20)
Now A = αV T . Each line of A is a sparse linear combination of V T 's lines, and the "active" atoms are optimally selected according to the data. The choice of the parameters (e i , a i ) 1≤i≤r and (η l ) 1≤i≤r is discussed in a forthcoming section.

C.2.4 A connection with graphs theory

In case a = 1, the matrix P e,a is the laplacian of an undirected fully connected and weighted graph with p nodes 1 . . . p, such that the weight of the vertex connecting a node i to a node j is

1 u i -u j e
2 [START_REF] Anderson | [END_REF]. As proposed in spectral graph theory [START_REF][END_REF]], this gives a natural interpretation of P e,a (and Q e,a ) eigenvectors as harmonic atoms in the graph's geometry. Each line of the matrix A can be seen as a function defined on a family of graphs determined by the observations locations, so that we enforce the regularity of A's lines according to the graphs geometry. Our approach is thereby close to the spectral graphs wavelets framework [START_REF] Hammond | [END_REF]]. However, the graphs wavelets are built on a single graph and a scaling parameter allows one to derive wavelets atoms corresponding to spectral bands of different sizes. In our case, the scales diversity is accounted for by building a dictionary of harmonics corresponding to different graphs. Indeed, as e increases, the weight associated to the most distant nodes (in the sense of u i -u j 2 ) becomes negligible, which implies that the corresponding graph laplacian is determined by nearby nodes, yielding "higher" frequencies harmonics.

C.3 The smoothness constraint on S

As previously mentioned, each PSF is a structured image. We can account for this through a sparsity constraint which this has proven effective in multiple frame PSFs super-resolution in the previous chapter.

Since we do not estimate individual PSFs directly, we instead constraint the eigen PSFs which are S's columns. Specifically, we promote S's columns sparsity with respect to a chosen dictionary Φ s . By definition, a typical imaging system's PSF concentrates most of its power in few pixels. Therefore a straightforward choice for Φ s is I n . In other words, we will enforce the sparsity of S's columns in the pixels domain.

On the other hand, we take Φ s as the second generation Starlet forward transform [Starck et al. 2011], without the coarse scale. The power of sparse prior in wavelet domain for inverse problems being well established, we shall online emphasize the fact that this particular choice of wavelet is particularly suitable for images with nearly isotropic features.

C.4 Algorithm

We define the sets

Ω 1 = {α ∈ M r,N (R)/ α[l, :] 0 ≤ η l , l = 1 . . . r} and Ω 2 = {(S, α) ∈ M nr (R) × M r,N (R)/SαV T ≥ Mnp(R) 0}.
The aforementioned constraints leads us to the following optimization problem :

min α,S 1 2 Y -F(SαV T ) 2 F + r i=1 w i Φ s s i 1 + ι Ω 1 (α) + ι Ω 2 (S, α). (V.21)
where denotes the Hadamard product and ι C denotes the indicator function of a set C (see B). The 1 term promotes the sparsity of S columns with respect to Φ s . The vectors (w i ) i weight the sparsity against the other constraints and allow some adaptivity of the penalty, with respect to the uncertainties propagated to each entry of S (see Chapter IV).

The parametric aspects of this method are made clear in the subsequent sections.

The Problem V.21 is globally non-convex because of the coupling between S and α and the l 0 constraint. In particular, the feasible set

{(S, α) ∈ M nr (R) × M r,N (R)/SαV T ≥ 0}, with N = rp is non-convex.
Therefore, one can at most expect to find a local minimum. To do so, we consider the following alternating minimization scheme :

1. Initialization : α 0 ∈ Ω 1 , with N = rp, S 0 = argmin S 1 2 Y - F(Sα 0 V T ) 2 F + r i=1 w i Φ s S[:, i] 1 s.t. Sα 0 V T ≥ 0 2. For k = 0 . . .k max : (a) α k+1 = argmin α 1 2 Y -F(S k αV T ) 2 F s.t. α[l, :] 0 ≤ η l , l = 1 . . . r, (b) S k+1 = argmin S 1 2 Y -F(Sα k+1 V T ) 2 F + r i=1 w i Φ s S[: , i] 1 s.t. Sα k+1 V T ≥ 0 .
The problem (a) remains non-convex ; yet there exists heuristic methods allowing one to approach a local minimum [Soussen et al. 2015, Blumensath & Davies 2008, Cartis & Thompson 2015]. The problem (b) is convex and can be solved efficiently. One can note that there is no positivity constraint in the sub-problem (a). This choice is motivated by two facts :

-the feasible set of (b) is non-empty for any choice of α k+1 ; -allowing α to be outside of the global problem feasible set (for S fixed) brings some robustness regarding local degenerated solutions.

There is an important body of work in the literature on alternate minimization schemes convergence, and in particular in the non-convex and non-smooth setting (see [START_REF] Bolte | [END_REF]] and the references therein). In the proposed scheme, the analysis is complicated by the asymmetry of the problems (a) and (b).

We define the function

H(α, S) = 1 2 Y -F(SαV T ) 2 F + r i=1 w i Φ s s i 1 (V.22)
and the matrix

S k = argmin S 1 2 S -S k 2 2 s.t. Sα k+1 V T ≥ 0. One immediate
sufficient condition for the sequence (H(α k , S k )) k to be decreasing (and thereby convergent) is

H(α k+1 , S k ) ≤ H(α k , S k ) (V.23)
which occurs if (S k , α k+1 ) stays sufficiently close to Ω 2 . Although we do not prove this always holds true, we observe on examples that the matrix S k α k+1 V T in general only has a few and small negative entries for k ≥ 1. This follows from the adequacy of the dictionary V for sparsely describing A's lines.

The complete method is given in Algorithm V.1. The resolution of the minimization sub-problems is detailed in Appendix D.

Algorithme V.1 Resolved components analysis (RCA)

1: Parameters estimation and initialization : Harmonic constraint parameters (e i , a i ) 1≤i≤r → V, A 0 Noise level, A 0 → W 0,0 2: Alternate minimization 3: Pour k = 0 to k max faire 4:

Pour j = 0 to j max faire 5:

S k = argmin S 1 2 Y -F(SA k ) 2 F + r i=1 W k,j [:, i] Φ s S[: , i] 1 s.t. SA k ≥ 0 6: update : W k,0 , S k → update(W k,j+1 ) 7: α k+1 = argmin α 1 2 Y -F(S k αV T ) 2 F s.t. α[l, :] 0 ≤ η l 8:
update : Noise level, α k+1 → W k+1,0 9:

A k+1 = α k+1 V T 10: A k+1 [i, :] = A k+1 [i, :]/ A k+1 [i, :] 2 , for i = 1 . . . r 11: Return : S kmax , A kmax .

C.5 Parameters setting

C.5.1 Components sparsity parameters

We consider the terms of the form w k,j Φ s s 1 , where k is the alternate minimization index and j is the re-weighted 1 minimization index. We first suppose that Φ s = I n . We decompose w k,j as :

w k,j = κβ k,j λ k (V.24)
Let consider the minimization problems in S in Algorithm V.1. Assuming that we simply minimize the quadratic term using the following steepest descent update rule,

S m+1 = S m + µF * (Y -F(S m A k ))A T k , (V.25)
for a well chosen step size µ, F * being the adjoint operator one can estimate the entry-wise standard deviations of the noise which propagates from the observations to the current solution S m+1 . For a given matrix X in M np (R), we assume that F takes the following general form

F(X) = [M 1 X[:, 1], . . . , M p X[: , p]].
We define F 2 as :

F 2 (X) = [(M 1 M 1 )X[:, 1], . . . , (M p M p )X[:, p]] (V.26)
We note B the observational noise (or model uncertainty) that we assume to gaussian, white and centered. The propagated noise entry-wise standard deviations are given by

Σ k = µ F 2 * (Var(B))(A T k A T k ), (V.27)
where Var() returns entry-wise variances and F 2 * is the adjoint operator of F 2 . Now one can proceed to a hypothesis testing on the signal presence in each entry of S m+1 based on Σ k [START_REF] Starck | [END_REF], and denoise S m+1 accordingly.

For instance, we define the noise-free version of S m+1 as follows :

Ŝm+1 [i 1 , i 2 ] = 0, if |S m+1 [i 1 , i 2 ]| ≤ κΣ k [i 1 , i 2 ] S m+1 [i 1 ,i 2 ] |S m+1 [i 1 ,i 2 ]| (|S m+1 [i 1 , i 2 ]| -κΣ k [i 1 , i 2 ]). otherwise; (V.28)
where κ controls the false detection probability ; the noise being gaussian, we typically choose 3 or 4 for κ.

The sequence ( Ŝm ) converges to a solution of the problem

argmin S 1 2 Y -F(Sα k U T ) 2 F + r i=1 κ λ k [:, i] S[:, i] 1 , (V.29) for λ k = κ/µΣ k .
One can find some material on minimization schemes in Appendix D. This choice yields a noise-free but biased solution because of the thresholding ; this is a well-known drawback of 1 norm based regularizations. The purpose of the vector β k,j is to mitigate this bias [Candès et al. 2008]. β k,0 is a vector with ones at all entries. At the step 6 in Algo V.1, β k,j is calculated as follows :

β k,j = 1 1 + |S k | κλ k , (V.30)
where all the operations are entry-wise and |S k | is the vector made of elementwise absolute values of S k entries. Qualitatively, this removes the strongest features from the 1 norm terms by giving them small weights, which makes the debiasing possible ; conversely, the entries dominated by noise get weights close to 1, so that the penalty remains unchanged.

For Φ s = I n we follow the same rational. To set the sparsity in the transform domain according to the noise induced uncertainty, we need to further propagate it (the noise) through the operator Φ s . Formally, we need to estimate the element-wise standard deviations of µΦ s F * (B)A T k . Let consider the intermediate random matrix Y F = F * (B). Assuming that F(F * (.)) = λId(.), (V.31) Y F 's lines are statistically independent. Therefore, within a given column of Y F A T k , the entries are statistically independent from one another. We deduce that the element-wise standard deviations of µΦ s F * (B)A T k are given by

Σ k = µ (Φ s Φ s )F 2 * (Var(B))(A T k A T k ). (V.32)
Then λ k is obtained as previously and β k,j is calculated as

β k,j = 1 1 + |ΦsS k | κλ k . (V.33)
The property V.31 is approximately true in the case of super-resolution.

C.5.2 Number of components

We do not propose a method to choose the number of components r. Yet, we observe that because of the sparsity constraint, some lines of the matrix α k+1 at the step 8 in Algorithm V.1 are equal to the null vector, when the number of components is overestimated. The corresponding lines in A k+1 and subsequently the corresponding columns in S k are simply discarded. This provides an intrinsic mean to select the number of components. Thus in practice, one can choose the initial r as the data set dimensionality from the embedding space point of view, which can be estimated based on a principal component analysis.

C.5.3 Proximity constraint parameters

In this section, we consider the functionals Ψ e i ,a i and especially the choice of the parameters e i and a i . Let assume that we have determine a suitable range for the parameters : (e i , a i ) ∈ S = [e min , e max ]×[a min , a max ] for i = 1 . . . r.

For a particular (e, a) we consider the matrix Q e,a and its eigenvectors matrix V e,a introduced in section C.2.3. As previously stated, we want the weights matrix A lines to be sparse with respect to Q e,a 's eigenvectors. In order to choose the parameters and initialize the weights matrix, we use the following greedy procedure. We consider a sequence of matrices (R i ) 1≤i≤r , with R 1 = Y. For i ∈ 1, r we define

J e,a (R i ) = max k∈ 1,p R i V e,a [:, k] 2 , (V.34)
and we note v * e,a the optimal eigenvector. We choose the i th couple of parameters as :

(e i , ai) = argmax (e,a)∈S J e,a (R i ). (V.35) A 0 [i, :] = v * e i ,a i and R i+1 = R i -R i V e,a
V T e,a . Regarding the set S, we choose the interval a min = 0 and a max = 2. This range allows the notch structure, assuming that e min ≥ 0 ; for a < 0, h e,a behaves as a low pass filter. For a ≥ 0, we observe that h e,a becomes a notch filter, with a notch frequency close to the null frequency for a ≥ 2. As previously stated, e determines the influence of two samples on one another corresponding coefficients in the matrix A in the algorithmic process. According to Section C.2.2, we set e min = 1. Let consider the graph G e introduced in section C.2.4. The higher is e, the lower is G e connexity. Considering that we are looking for global features (yet localized in the field frequency domain), the highest possible value of e should guarantee that the graph G e is connected. This gives us a practical upper bound for e. Once S is determined, we discretize this set, with a logarithmic step, in such a way to have more samples close to (e min , a min ) which correspond to low notch frequencies. We solve approximately Problem V.35 by taking the best couple of parameters in the discretized version of S. This step is the most computationally demanding, especially for large data samples.

C.5.4 Weights matrix sparsity parameters

The parameters η l are implicitly set by the minimization scheme used at step 8 in V.1. This is detailed in Appendix D.

D Numerical experiments

In this section, we present the data used to test the proposed method, the simulation realized and comparisons to other existing methods for both dimensionality reduction and super-resolution aspects.

D.1 Data

The data set consists of simulated optical Euclid PSFs as in Chapter IV, for a wavelength of 600µm. The PSFs distribution across the field is shown on Fig.

V.3. These PSFs account for mirrors polishing imperfections, manufacturing and alignment errors and thermal stability of the telescope.

D.2 Simulation

We applied different dimension reduction algorithms to a set of 500 PSFs located in the blue box on Fig. V.3. We applied the algorithms to different observations of the fields, with varying level of white gaussian noise. For a discrete signal s of length N corrupted with a white gaussian noise b, we define the signal to noise ratio (SNR) as :

SNR = s 2 2 N σ 2 b . (V.36)

D.3 Quality assessment

In astronomical surveys, the estimated PSF's shape is particularly important ; precisely, one has to be able to capture the PSF anisotropy. We recall that for an image X = (x ij ) i,j , the central moments are defined as

µ p,q (X) = i j (i -i c ) p (j -j c ) q x ij (V.37)
with (p, q) ∈ N 2 , (i c , j c ) are the image centroid coordinates. The moments µ 2,0 and µ 0,2 quantifies the light intensity spreading relatively to the lines {(i c , y), y ∈ R} and {(x, j c ), x ∈ R} respectively. Now we consider the moment µ 1,1 . We introduce the centered and rotated pixels coordinates (x i,θ , y j,θ ) defined by the system of equations

x i,θ cos(θ) + y j,θ sin(θ) = i -i c (V.38) -x i,θ sin(θ) + y j,θ cos(θ) = j -j c , (V.39)
for some θ ∈ [0, 2π]. Then we have

µ 1,1 = i j [ sin(2θ) 2 (-x 2 i,θ + y 2 j,θ ) + (2 cos 2 (θ) -1)x i,θ y j,θ ]x ij , (V.40)
and in particular,

µ 1,1 = i j [ 1 2 (-x 2 i, π 4 + y 2 j, π 4
)]x ij . It becomes clear that µ 1,1 quantifies the light intensity spreading with respect to the pixels grid diagonals.

The ellipticity parameters are defined as,

e 1 (X) = µ 2,0 (X) -µ 0,2 (X) µ 2,0 (X) + µ 0,2 (X) (V.41) e 2 (X) = 2µ 1,1 (X) µ 2,0 (X) + µ 0,2 (X) . (V.42)
We define the vector γ(X) = [e 1 (X), e 2 (X)] T . This vector characterizes how much X departs from an isotropic shape and indicates its main direction of orientation. It plays a central theoretical and practical role in weak lensing based dark matter characterization [START_REF] Dodelson | Modern cosmology[END_REF]].

Another important geometric feature is the so-called PSF size. It has been shown that the size error is a major contributor to the systematics in weak gravitational lensing surveys [START_REF] Paulin-Henriksson | [END_REF]. We characterize the size of a PSF X as follows :

S(X) = ( i j ((i -i c ) 2 + (j -j c ) 2 )x ij i j x ij ) 1/2 . (V.43)
Assuming that a given PSF is a 2D discrete probability distribution, this quantity measures how much this distribution is spread around its mean [i c , j c ] T . Let note (X i ) 1≤i≤p the set of "original" PSFs and ( Xi ) 1≤i≤p the set of corresponding estimated PSFs with one of the compared methods, at a given SNR. The reconstruction quality is accessed through the following quantities :

-the average error on the ellipticity vector :

E γ = p i=1 γ(X i ) - γ( Xi ) 2 /p ; -noting Γ = [γ(X 1 ) -γ( X1 ), . . . , γ(X p ) -γ( Xp )],
the dispersion of the errors on the ellipticity vector is measured through the nuclear norm B γ = Γ * ;

-the average absolute error on the size :

E S = p i=1 |S(X i ) -S( Xi )|/p in pixels ;
-the dispersion of the errors on the size : σ S = std((S(X i ) -S( Xi )) i ), in pixels.

E Numerical experiments

E.1 Results

E.1.1 Dimension reduction

In this section, we compare RCA to PCA, GMRA and the software PS-FEx. We ran a PCA with different number of principal components between 0 and 15. 10 was the value which provided the best results. GMRA input was the data set intrinsic dimension [Little 2011], two, since the PSFs only vary as a function of their position in the field ; we provided the absolute squared quadratic error allowed with respect to the observed data based on the observation noise level. For PSFEx, we used 15 components. Finally, RCA used up to 15 components, and effectively, 2 and 4 components respectively for the lowest SNR fields realization. As previously mentioned, we assess the components sparsity's constraint :

-on the one hand we consider Φ s = I n which enforces the components sparsity in pixels domain ; this is referred to as "RCA" in the plots ;

-on the other hand, we take Φ s as the second generation Starlet forward transform [Starck et al. 2011], without the coarse scale ; this is referred to as "RCA analysis" in the plots.

One can see on the left plot in Fig.

V.4 that the proposed method is at least 10 times more accurate on the ellipticity vector than the other considered methods. Moreover the right plot shows that the accuracy is way more stable. This is true for both choice of the dictionary Φ s . Fig.

V.9 shows that the estimated size S( Xi ) is very sensitive to the choice of the dictionary Φ s . The results are by far more accurate with a sparsity constraint on the components in wavelet domain than in direct domain.

For a given estimate of the PSF at a given location, the error on the size parameter is more sensitive to errors on the core of the PSF (main lobe and first rings) and less sensitive to errors on the outer part of the PSF than one would expect regarding the error on the ellipticity vector. The error on the outer part of the PSF is essentially related to the observational noise, whereas the error on core of the PSF -which has a high SNR -is more related to the method induced bias. This explains why the PCA performs quite well for this parameter. On the other hand, the bias induced by the sparsity is not only related to the dictionary choice, but also to the underlying data model with respect to the chosen dictionary. As previously explained, the components sparsity term is set in such a way to penalize any feature which does not emerge from the propagated noise, which is a source of bias. By using wavelets, we might recover features which are dominated by noise in pixel domain as long as the wavelet filters profile at given scale and direction, matches those features spatial structure. Thus, we expect less error on the reconstructed PSF's core by using wavelets.

We might also consider two distinct ways of using sparsity for the components :

-we can model each component as s = Φ T s α, with α sparse, which is known in the sparse recovery literature as synthesis prior ;

-we can alternately constraint Φ s s to be sparse. This priors are equivalent if the dictionary is unitary [Elad et al. 2007]. Therefore the pixel domain sparsity constraint can be considered as falling into both framework. However, the two priors are no longer equivalent and potentially yields quite different solutions for overcomplete dictionaries.

We observe in practice that unless the simulated PSFs are strictly sparse with respect to the chosen dictionary -this includes redundant wavelet dictionaries, the synthesis prior yields a bias on the reconstructed PSF size, since the estimated PSFs are sparse linear combinations of atoms which are in general sharper than a typical PSF profile. The analysis prior is somehow weaker and appears to be more suitable for approximately sparse data.

We do not observe a significant difference between these methods with respect to the mean squared error, except for GMRA which gave noisier reconstructions.

We applied the aforementioned methods to the PSFs field previously used, with additional 30 corners PSFs and 30 localized PSFs as shown on Fig. One can observe that the proposed method gives noiseless and rather accurate PSFs reconstruction for both the center, the corners and the localized area of the field (see Fig. V.3). One can also see that we fail to capture accurately the rings pattern in the corners and the localized area. The dictionary Φ s considered are not specifically adapted to curve-like structures. The ring patterns varies across the FOV but are locally correlated. Therefore, they can only be recovered where the PSFs are sufficiently dense and numerous, which is the case at the FOV's center.

PCA and PSFEx yield a significant increase of the SNR in their estimated PSFs at the center and in the localized area. Yet, they fail to do so in the Finally, the poor results obtained with GMRA can be explained by the fact that the underlying manifold sampling is not sufficiently dense for the tangent spaces to be estimated reliably.

E.1.2 Super-resolution

In this section, the data are additionally downsampled to Euclid telescope resolution. PCA and GMRA does not handle the downsampling. Therefore we only consider PSFEx and RCA in this section. For each method, we estimate an upsampled version of each PSF, with a factor 2 in lines and columns ; in case of Euclid, this is enough to have a Nyquist frequency greater than half the signal spatial bandwidth [Cropper 2013].

As previously, RCA Analysis refers to the proposed method, with the dic- tionary Φ s chosen as the second generation Starlet forward transform [Starck et al. 2011], without the coarse scale ; RCA LSQ refers to the proposed method with the dictionary Φ s chosen as the identity matrix, and the weight matrix A simply calculated as

A = argmin A 1 2 Y -F( SA) 2 F , (V.44)
S being the current estimate of the components matrix. Among all the methods previously considered for comparison, PSFEx is the only one handling the undersampling. As for the dimension reduction experiment, the proposed method with Φ s chosen as a wavelet dictionary is at least one order of magnitude more accurate over the shape parameters and the mean square error. Besides, Fig.

V.10 shows that the proximity constraint over the matrix A allows one to select a significantly better optimum than a simple least square update of A. Indeed, regularizing the weight matrix estimation reinforces the rejection of F's null space.

As previously, we restored the complete field of Fig.

V.3 for a linear SNR of 40, using "RCA Analysis", with undersampled input PSFs as shown in Fig.

V.11.

The figure V.12 shows consistent results with the dimension reduction experiment. In particular, the corners PSFs restoration is obviously more accurate.

F Reproducible research

In the spirit of participating in reproducible research, the data and the codes used to generate the plots presented in this paper will be made available at http://www.cosmostat.org/software/rca/.

G Conclusion

We introduced RCA which is a dimension reduction method for continuous and positive data field which is noise robust and handles undersampled data. As a linear dimension reduction method, RCA computes the input data as linear combinations of few components which are estimated, as well as the linear combination coefficients, through a matrix factorization. The method was tested over a field of simulated Euclid telescope PSFs. We show that constraining both the components matrix and the coefficients matrix using sparsity yield at least one order of magnitude more accurate PSFs restoration than existing methods, with respect to the PSFs shapes parameters. In particular, we show that the analysis formulation of the sparsity constraint over the components is particularly suitable for capturing accurately the PSFs sizes. We also show that constraining the coefficients matrix yields a significantly better identification of the PSFs embedding subspace when the data are undersampled. Chapitre VI This chapter focuses on interpolating a PSF field at arbitrary locations from a set of perfectly known PSFs spread out across the fov. Assuming that the PSF can be modeled locally as a convolution kernel, we use "PSF field" to refer to a continuous function mapping each point in the instrument focal plane surface to the corresponding convolution kernel. This problem has driven a lot of attention within the astronomers community over the last decade [Gentile et al. 2013], especially due to the strengthening of accuracy requirements in the recent and future spatial surveys. We review some of the state-of-the-art method in the next section.

PSF field interpolation

A PSF field interpolation state-of-the-art

Most of the PSF interpolation methods fall into one of the two following categories :

• optical model based methods which derives the PSF from a physical model of the instrument itself ;

• data driven methods which estimate the PSFs at given locations in the fov based on local measurements extracted from real images.

In the first category, we can mention the work in [START_REF] Jarvis | [END_REF]. The authors built a parametric model of the PSF shape and size by modeling individually the main physical causes which make the PSF variable ; they found that this model is able to reproduce a substantial part of the PSF anisotropy and size for the ground based telescope Blanco telescope. In [Stabenau et al. 2007], the authors simulate the PSF time variations due to thermal drift, jitter and structural vibrations, using the method of ray tracing, for the SNAP telescope. However, these methods require making assumptions on the instrument physic, and not all the phenomena at play in the image forming process can be satisfactorily modeled. Thus, we only consider model-free methods for comparisons in the numerical experiments.

The data driven approaches rely on the fact that in practice, unresolved stars images in the observed field give a measurement of the PSFs at those stars locations. One recurrent scheme consists in first expanding the observed unresolved stars over some given analytic function basis. This gives a more compact representation of the data (and potentially allows denoising). The PSF field is then obtained through a 2D polynomial fitting typically, of the representation coefficients for each element of the chosen basis. In [Romano et al. 2010], this scheme is applied to model the PSF of the Large Binocular Cameras, using the Shapelets basis ( [Refregier 2003, Massey & Refregier 2005]). However, it might not be possible to capture all the PSF structures through a finite expansion over an analytic function basis. One can instead learn the representation basis from the data themselves. For example in [START_REF] Jee | [END_REF]], a principal component analysis (PCA) is performed over hundreds archival images of stars from the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS), from different exposures. The optimal representation basis is then chosen as the first principal components and the representation coefficients are fitted in each exposure separately with a bivariate polynomial. In [Jarvis & Jain 2004], the PSF modeling also relies on a PCA, not of the images themselves, but instead of a set of features (for example a shape parameter). In the two aforementioned methods, fitting a global polynomial model can yield an oversmoothing of the PSF field, or instead a Runge phenomenon, if the polynomial's degree is too high [Epperson 1987]. In [Gentile et al. 2013], the PSF modeling is recast as a spatial interpolation problem. In general statistics, the spatial interpolation methods aim at interpolating a field exploiting the spatial correlation of the data. Some of these methods, for example the Kriging, are presented and applied to PSF interpolation in the previously mentioned work. We consider the two best performing approaches tested in [Gentile et al. 2013] for comparisons in the numerical experiments. More recently a PSF interpolation scheme has been proposed in [START_REF] Suksmono | Reconstruction of Complex-Valued Fractional Brownian Motion Fields Based on Compressive Sampling and Its Application to PSF Interpolation in Weak Lensing Survey[END_REF]] based on, compressive sampling (CS) ideas. The method aims at computing the PSF ellipticity paramters at every points over an uniform grid using the available randomly distributed measurements. This yields an ill-posed inverse problem which is regularized by assuming that the shape parameters constitute a compressible field in frequency domain. However, this approach has a narrow scope of applications since the "full" PSF itself is not estimated.

Let finally mention a somehow hybrid approach. Indeed in [Piotrowski et al. 2013], the authors propose to build an accurate polynomial model of the PSF of the "Pi of the Sky" telescope from detailed laboratory measurements of optical PSF, pixels sensitivity and pixels response. This model is then tuned to real sky data.

In this paper, we propose a non-parametric interpolation method that relies on geometric tools. We propose an intuitive framework for learning a PSFs set underlying geometry using Optimal Transport distances. In Section C, we introduce and motivate the mathematical framework ; we describe the proposed method and algorithm in Section D ; we present some numerical experiments and discuss the results in Section E ; the interpolation algorithm parameters and some numerical considerations as discussed in Section F.

B Notations

Each PSF is treated either as a matrix X = (x ij ) 1≤i≤n 1≤j≤n or as a vector in R n 2 , x = (x k ) 1≤k≤n 2 . For N ∈ N, we denote Σ N the set of permutations of 1, N . We denote P(R d ) the set of discrete probability measures over R d . We denote 1 n the column vector of n ones. We denote I p the identity matrix of size p × p.

C Numerical Optimal Transport

C.1 Motivation

We consider a field of monochromatic PSFs. The only factor of variability is the unresolved object's position in the fov. One can see from Fig. VI.1 that PSFs shapes vary non-linearly across the field. This is related to differences in the paths an electromagnetic wave would follow from the telescope pupil entrance to the focal plane, depending on the source position ; this results in different diffraction's figures, which are the observed PSFs. This can be interpreted as the curvature in R n 2 of the PSFs manifold, occurring between distant PSFs in the fov. As Fig. VI.2 shows, we want to combine available PSFs measurements given by unresolved objects images to estimate the PSFs at galaxies locations in the fov.

Optimal Transport (OT) appears to be suitable to tackle this problem. Indeed, it gives a way of measuring distances between PSFs which accounts for the aforementioned curvature.

Besides, as already highlighted in Chapter III, Section B, it comes with a recipe for computing geodesics which is interesting for interpolation. This is detailed in the next section. Moreover, it is classical to assume that a PSF has an unitary l 1 norm. This implies a constraint of mass conservation in the interpolation which is naturally integrated in the OT framework. In the next section, we precise the OT tools and concepts we use.

C.2 General notions

C.2.1 Sliced Wasserstein Distance

We consider two points matrices X = [x 1 , . . . , x N ] and Y = [y 1 , . . . , y N ] in R d×N , for some integers d and N . By abuse of language, we do not distinguish these matrices and the point clouds in R d they represent, where there is no risk of confusion. We define the distributions µ

X = N i=1 1 N δ x i and µ Y = N i=1
1 N δ y i , where for x and y in R d , δ x (y) = 1 if x = y and δ x (y) = 0 otherwise. The optimal assignment problem between µ X and µ Y consists in finding a permutation σ ∈ Σ N which minimizes

N i=1 x i -y σ(i) p , (VI.1)
where . is a norm in R d and p ≥ 1. Problem VI.1 can be rewritten as a particular instance of the Monge-Kantorovich mass transportation problem [START_REF] Kantorovitch | ON THE TRANSLOCATION OF MASSES[END_REF]. We note σ * a minimizer of Problem VI.1. It has been shown

that W p (µ X , µ Y ) = ( N i=1 x i -y σ * (i) p p )
1 p defines a distance on P(R d ) which is the so-called Wasserstein Distance p of the two distributions. The norm . p is often referred to as the ground metric and is an important parameter of the Transport Distances in practical applications. Examples of theoretical and computational applications of optimal transport can be found in [Papadakis et al. 2013].

Let explicit the notion of displacement interpolation in this discrete setting. Let consider a point x i in the first point cloud and the assigned point y σ * (i) in the second point cloud. Let γ i : [0, 1] → R d a curve verifying γ i (0) = x i and γ i (1) = y σ * (i) . The action on a material particle moving from x i to y σ * (i) along γ i between time s and time t might be defined as

A s,t p (γ i ) = t s γi (τ ) p dτ, (VI.2)
where γi (τ ) is the particle velocity at time τ . We define the set

Γ i = {γ : [0, 1] → R d , γ(0) = x i , γ(1) = y σ * (i)
}. An action minimizing path between x i and y σ * (i) is a solution of min

γ∈Γ i A s,t p (γ i ). (VI.3)
Interestingly, it has been shown that a geodesic path between the distributions µ X and µ Y in P(R d ) equipped with the Wassertein distance can be computed by advecting particles along action minimizing paths between pairs of assigned points in the two point clouds [Villani 2009]. Concretely, let γ i * denotes an action minimizing path between x i and y σ * (i) . For t ∈ [0, 1], we define the discrete distribution

µ t = N i=1 1 N δ γ i * (t) . (VI.4)
Then the parametric curve

Γ XY : [0, 1] → P(R d ), t → µ t (VI.5)
is a geodesic path in P(R d ) between the distributions µ X and µ Y . This realizes a displacement interpolation of the two considered distributions ; this powerful notion had first been introduced in [McCann 1997]. For a given t, the distribution µ t can be seen as the barycenter of the two end distributions with the barycentric weights 1 -t and t respectively in the Wasserstein metric. For d = 1, i.e. in the 1D case, σ * is known in closed form : if we consider two permutations σ X and σ Y in Σ N verifying

x σ X (1) ≤ • • • ≤ x σ X (N ) (VI.6) x σ Y (1) ≤ • • • ≤ x σ Y (N ) , (VI.7) then σ * = σ Y • σ -1 X , where σ -1 X is a permutation verifying σ -1 X • σ X = Id. For d > 1,
there is no known closed-form expression for σ * . Various methods have been proposed for solving Problem VI.1 ; in particular, it can be recasted as a linear program. Yet, the fastest known algorithms have a running time of O(N 2.5 log(N )) [Burkard et al. 2009], which is prohibitive for large scale image processing applications. This has motivated the introduction of the sliced Wasserstein Distance [Rabin et al. 2011, Bonneel et al. 2015], which consists in the sum of 1D Wassertein distances of the projected point clouds :

SW p (µ X , µ Y ) p = S d-1 W p (µ Xu , µ Yu ) p du, (VI.8) where p ≥ 1, S d-1 = {x ∈ R d / x = 1}, X u = {u T x i , i = 1...N } ⊂ R N and
Y u is similarly defined.

C.2.2 Sliced Wassertein Barycenter

We consider a family {X 1 , . . . , X K } of point clouds in R d : for k ∈ 1, K , X k = [x k1 , . . . , x kN ] ∈ R d×N . The associated distributions µ X k are defined as previously. A barycenter of these distributions in the Wasserstein metric is defined as

µ Bar = argmin µ K k=1 w k W p (µ, µ X k ) p , (VI.9)
for some positive weights w k verifying K k=1 w k = 1. Somehow, this generalizes the concept of displacement interpolation for an arbitrary finite number of distributions. In various applications Wasserstein barycenters appear to be more suitable than Euclidean or more sophisticated barycenters ( [Cuturi & Doucet 2013], [Gramfort et al. 2015]) ; in particular they are robust to shifts and "elastic" deformations. As previously, in 1D (i.e. d = 1) and for p = 2, this barycenter is known in closed-form. It is the uniform discrete probability measure defined over the set

{x Bari = K k=1 w k x k [σ k (i)], i = 1 . . . N }, (VI.10)
where σ k is a permutation in Σ N which sorts x k entries (i.e.

x k [σ k (1)] ≤ • • • ≤ x k [σ k (N )]).
For d > 1, computing this barycenter is in general intractable [START_REF] Gangbo | [END_REF]]. Therefore in this paper we consider the Sliced Wasserstein Barycenter defined by replacing the Wasserstein distance by its sliced counterpart in the Wasserstein barycenter definition :

µ SBar = argmin µ K k=1 w k SW p (µ, µ X k ) p .
(VI.11)

C.3 Data representation in the Transport Framework

We consider a N l × N c gray levels image. We denote x = (x i ) 1≤i≤N the vector of pixels intensities rearranged in lines lexicographic order, with N = N l N c . The image is treated as a point cloud defined as (VI.12) This point cloud can be viewed as a representation of the considered image as a discrete surface in R 3 . Given a set of unresolved objects images, the idea in the following section would be to match and or compute barycenters of the associated surfaces based on sliced transport. We use a ground metric C taking the following form : for two vectors p = [p 1 , p 2 , p 3 ] T and q = [q 1 , q 2 , q 3 ] T in R 3 ,

X = [v 1 , • • • , v N ], with v i = [x i , i -1 N c , (i -1(modN c ))] T .
C 2 (x, y) = (p 1 -q 1 ) 2 + β 2 * ((p 2 -q 2 ) 2 + (p 3 -q 3 ) 2 ), (VI.13)
where β is a strictly positive real. As a reminder, in the applications, the first components p 1 and q 1 will be pixels intensities and the two other components will be pixels positions. For two images x = (x i ) 1≤i≤N and y = (y i ) 1≤i≤N , let X and Y be the associated point clouds and σ * an optimal assignment of Y's points to X's points : then

W 2 (µ X , µ Y ) = N i=1 (x i -y σ * (i) ) 2 +β 2 * (( i -1 N c - σ * (i) -1 N c ) 2 + ((i -1(modN c )) -(σ * (i) -1(modN c ))) 2 ).
(VI.14) Thus, if σ * is the identity, W 2 (µ X , µ Y ) = xy 2 . Therefore, we can consider W 2 (and SW 2 ) as a generalization of the euclidian distance in the pixels domain.

Conversely, from a given point cloud Z ∈ R 3×N , we can go backward to an image by discretizing the 2D function Φ[Z] defined as

Φ[Z](x, y) = N i=1 Z[1, i]δ Z[2,i],Z[3,i] (x, y).
(VI.15)

D Transport based PSF field interpolation

We consider a set of K PSFs S = (x k ) 1≤k≤K , located at the positions (u k ) 1≤k≤K in the fov. We note (X k ) 1≤k≤K the associated point clouds as defined in Section C.3. We want to estimate the PSF at a new location u. We note this PSF x u and we note X u the associated point cloud :

x u ≡ Φ[X u ],
(VI.16)

where the operator Φ is defined in Eq.VI.15. We want to estimate X u as a Wassertein barycenter of a subset of PSFs indexed by I(u) in the neighborhood of u in the fov :

X u = argmin X k∈I(u) w k (u)SW 2 (µ X , µ X k ) 2 , (VI.17)
with w k (u) ≥ 0 and k∈I(u) w k (u) = 1. I(u) simply indexes the nearest neighbors of u in the set of locations (u k ) 1≤k≤K . The question of the size of this neighborhood is postponed to Section F.1.

The weights (w k (u)) k∈I(u) can be seen as generalized barycentric coordinates of X u relatively to the clouds (X k ) k∈I(u) , in the Wassertein metric.

These weights are calculated in three steps that we detail in the following sections :

-Embedding : we calculate the pairwise approximated Wasserstein 2 distances over the set (X k ) k∈I(u) and determine an euclidean embedding of (X k ) k∈I(u) that preserves the Wasserstein distances ; we get a set of coordinates (r k ) 1≤k≤|I(u)| , where |I(u)| is the number of elements in I(u) ;

-Interpolation : we estimate the representation r u of X u in the previously calculated embedding by interpolating the set (r k ) 1≤k≤|I(u)| coordinateswise ;

-Weights setting : the weights (w k (u)) k∈I(u) are calculated as the euclidean barycentric coordinates of r u relatively to (r k ) 1≤k≤|I(u)| .

D.1 Local non-linear dimension reduction

We recall that we want to estimate the PSF at a position u in the fov, given the PSFs located at the positions (u k ) 1≤k≤K . We consider the p nearest neighbors of u in (u k ) 1≤k≤K . This defines the set I(u) aforementioned.

As we will see in the Numerical experiments Section, the parameter β in Eq.VI.14 is set so that for the closest PSFs in the fov, the Wassertein distance is equal to the euclidian distance. Indeed, the euclidean distance is usually assumed to be a good approximation of the geodesic distance between close points on a given smooth Manifold (see for example [START_REF] Tenenbaum | [END_REF],Donoho & Grimes 2003]). However, the more two PSFs are distant in the fov, the more the PSFs manifold curvature manifests through the relative warping of their structures. By minimizing the amount of work needed to push one of the PSFs toward the other, the Wasserstein metric can keep track of the warping to a certain extent, thus unfolding locally the Manifold. For this reason, the Wasserstein distance can potentially give a faithful approximation of the geodesic distances, on broader neighborhoods than the euclidean distance. We define the p × p local pairwise distance matrix M W (u) as :

M W (u)[i, j] ≈ W 2 (µ X I(u)[i] , µ X I(u)[j] ), (i, j) ∈ 1, p 2 .
(VI.18)

Then we apply the Multidimensional Scaling (MDS) procedure to calculate the embedding [Abdi 2007]. The first step consists in converting the distances matrix M W (u) into a Gram matrix X i i.e. whose entries are given by X i [i, j] = r i , r j , so that M W (u)[i, j] = r i -r j 2 . The vector r i is the low dimensional embedding of X I(u) [i] .

For this purpose we introduce the "centering matrix" [Marden 1995] defined as

C = I p - 1 p 1 p 1 T p , (VI.19)
where 1 p is a row vector made of p ones. Multiplying a matrix by C on the left has the effect of subtracting each row its mean ; C being symmetric, multiplying a matrix by C on the right subtracts each lines its mean.

The cross-product matrix is calculated as (VI.20) where M W (u) is squared entry-wise. We define the embedding coordinates matrix R = [r 1 , . . . , r n ]. Assuming that such an embedding exists and R has null lines and columns means, one can show that X i = R T R (see appendix F). Under this hypothesis, X i is a symmetric matrix and therefore can be orthogonally diagonalized :

X i = - 1 2 CM W (u) 2 C,
X i = V T SV, (VI.21)
where V T V = VV T = I p and S is a non-negative diagonal matrix. It follows that R can be calculated as

R = QS 1 2 V, (VI.22)
where Q is an orthogonal matrix that we set to the identity in practice [START_REF] Dokmanic | [END_REF]. Depending on neighborhood's size p, the last diagonal values of S can be neglected so that the last lines of V can be discarded ; besides this represents a way of analyzing locally the Manifold dimensionality [START_REF] Tenenbaum | [END_REF]].

D.2 Field-of-view mapping

In this section, we estimate the local low dimensional embedding of the unknown PSF located at u in the fov. We note this embedding r u . We note d the dimension of the vectors (r i ) 1≤i≤p determined in the previous section, and r u . To compute the i th component of r u we determine an interpolating function

f i : R 2 → R /f i (u j ) = r j [i], ∀j ∈ 1, p . (VI.23)
This is a standard surface interpolation problem that we solve using the socalled thin-plate spline [Eberly 2002] which is appealing because of its physical interpretation : it is the exact interpolating function that minimizes the "bending energy" defined as :

E(f ) = R 2 ( ∂ 2 f (x, y) ∂ 2 x ) 2 + 2 * ( ∂ 2 f (x, y) ∂x∂y ) 2 + ( ∂ 2 f (x, y) ∂ 2 y ) 2 .
(VI.24) Thus, the function f i takes the following form :

f i (x) = p j=1 a ij x -u j 2 2 ln( x -u j 2 ), (VI.25)
and the coefficients a ij are calculated so that f i takes the prescribed values at the control points (u j ) 1≤j≤p . This way, each component of r u can be estimated.

D.3 Barycentric coordinates

The embedding of the PSF at the position u in the fov has been determined. Now, we want to estimate the point cloud X u from this embedding and the neighbor PSFs. In other terms, we need to determine the weights w k (u) in Eq. VI.17. In Section D.1, we calculated an isometric embedding of the matrices (X k ) k∈I(u) . Because of this isometry, for a set of positive weights (w i ) 1≤i≤p verifying p i=1 w i = 1, the problem (VI.27) This somehow goes along the same line as the Locally Linear Embedding [Roweis & Saul 2000] and it provides a mean for computing X u from r u . Indeed, we first can consider the following barycentric coordinates problem : We note (w * i ) 1≤i≤p the optimal tuple and we define r * u = p i=1 w * i r i . r u ≈ r * u and r * u is solution of the Problem VI.26, with the weights (w * i ) 1≤i≤p . Therefore, we compute X u as

w i SW 2 (µ X , µ X I(u)[i] ) 2 .
X u = argmin X p i=1 w * i SW 2 (µ X , µ X I(u)[i] ) 2 . (VI.29)
This approach is in a way similar to the procedure described in [START_REF] Gower | [END_REF]].

Besides, the task of estimating barycentric coordinates in the Wasserstein metric has been recently addressed in [Bonneel et al. 2016], in a different numerical Optimal Transport framework. The final image is obtained by discretizing the function Φ[X u ] defined by Eq.VI.15. This step is discussed in Section F.

D.4 Algorithm

We want to estimate D PSFs at random locations (v i ) 1≤i≤D from a set of K PSFs at known locations in the fov (see the Section D introduction). The whole procedure is summarized in Algorithm V I.1. Form the pairwise Wasserstein distances matrix over the set (X j ) j∈I(v i ) and compute a local euclidian embedding ; the results is a set of vectors (r j ) 1≤j≤p in R d

8:

Estimate the embedded coordinates at the location v i using a thin-plate spline interpolation coordinate-wise ; this results in a vector r v i 9:

Compute the barycentric coordinates of r v i relatively to the vectors (r j ) 1≤j≤p 10:

Compute the approximate Wassertein barycenter of (X j ) j∈I(v i ) using the previously calculated barycentric coordinates as weights 11:

Compute the interpolated PSF x v i from this barycenter (see Eq. VI.15) 12: Return :

(x v i ) 1≤i≤D .
This procedure is generic in the sense that the Wassertein metric can be replace by an arbitrary metric provided that one is able to compute geodesic distances and geodesics as we will see in Section F.

E Numerical results

We tested the proposed method on a set of 550 simulated Euclid telescope optical PSFs. The PSFs are distributed in the fov according to Fig. VI.3. We split the data into a "learning set" made of 300 observed PSFs and a "test set" corresponding to 250 unknown PSFs.

E.1 Quality assessment

As in the previous chapters, the quality assessment are the PSF's shape parameters. We recall their definitions thereafter. The central moments of an image X = (x ij ) i,j are defined as :

µ p,q (X) = i j (i -i c ) p (j -j c ) p x ij (VI.30)
with (p, q) ∈ N 2 , (i c , j c ) being the image centroid coordinates.

The ellipticity parameters are defined as follows :

e 1 (X) = µ 2,0 (X) -µ 0,2 (X) µ 2,0 (X) + µ 0,2 (X) (VI.31) e 2 (X) = 2µ 1,1 (X) µ 2,0 (X) + µ 0,2 (X)
.

(VI.32)

The vector γ(X) = [e 1 (X), e 2 (X)] T tells how much X departs from an isotropic distribution of luminosity and gives one the main direction of elongation.

The PSF's "size" is defined as

S(X) = ( i j ((i -i c ) 2 + (j -j c ) 2 )x ij i j x ij ) 1/2 . (VI.33)
We note (Im i ) 1≤i≤D the set of tests PSFs and ( Im i ) 1≤i≤D the set of corresponding PSFs interpolated with a given method. The reconstruction quality is accessed through the following quantities :

-the average error on the ellipticity vector :

E γ = D i=1 γ(Im i ) - γ( Im i

E.2 Experiments

We compare the proposed method to the Inverse Distance Weighting (IDW) and Radial Basis Function (RBF) based interpolation. These methods are described in [Gentile et al. 2013] as the one which performs the best on the considered PSFs interpolation. For the first method, the weights are calculated based on squared distances in the fov ; this is often referred as "Inverse distance-squared interpolator". We used the radial basis function f (r) = rln(r) ; thus the second method is nothing but a thin plate spline interpolation as described in Section D.2. The 300 observed PSFs are first decomposed using a principal components analysis and retaining 40 principal components. This yields a representation with 40 coefficients for each PSF. These coefficients are interpolated components-wise with the two comparison methods and the interpolated PSFs are derived. As the method proposed, the IDW and RBF based interpolation methods are local in the sense that a PSF is interpolated at given location in the fov based on a given number of PSFs observed in the vicinity of the considered location. Thus, the three methods will be applied using different numbers of "neighbors PSFs".

The results are shown in Fig. VI.4,VI.5,VI.6. Interestingly, we observed quite different behaviors with respect to the different criteria. This will be discussed in the next section. However, the less accurate method appears to be the IDW based method. Indeed, it gives the less accurate results in terms of shape ; in particular, it is almost an order of magnitude less accurate than the two other methods for the ellipticity. Furthermore, it is almost constantly less accurate than TraIn with respect to the pixels mean square error. The RBF based method and the proposed one gives comparable errors with respect to the shape parameters ; it is slightly more accurate on the ellipticity and slightly less accurate on the shape. The lowest errors with the RBF based method are obtained for a number neighbors around 15 for the shape parameters. The proposed method has a quite stable accuracy with respect to the number of neighbors for all the criteria. Moreover, it is typically several orders of magnitude more accurate than the RBF based one on the pixels mean square errors for a number of neighbors smaller than 15. Therefore, the TraIn method is globally the most accurate. We give some examples of absolute error images i.e.

E.3 Discussion

In this Section, we discuss the performances of the tested methods. We consider the q first principal component derived from the learning PSFs set. In the numerical experiments, we set q = 40 for the RBF based methods. The Figure VI.8 illustrates how coefficients magnitudes vary across the field for the whole data set (including the test PSFs) for some of the principal compo- nents. The two first components exhibit a smooth spatial evolution so that it is possible to interpolate accurately the corresponding surfaces with RBFs using few observations or control points in the vicinity of the interpolation point. This smoothness can be seen in the scatter plot associated with this two first components. Indeed, one needs a smooth warping to transform the scatter plot distribution into the PSFs spatial distribution. On the contrary, the magnitudes have sharper variations for the two other components. This can also be seen from the scatter plot, since a highly non-linear transform would be needed for mapping the scatter plot distribution into the PSFs spatial distribution.

In this case, more control points are required for the RBF interpolation to be robust to fast local variations. The scatter plots of Fig. VI.8 are nothing but orthogonal projections of the PSFs manifold over the vector plans spanned by the respective pairs of principal components. As the indexes of the principal components increase, the manifold complexity becomes more apparent and manifests in fast spatial variations of the representation coefficients. Thus, when a few control points are used, the RBF interpolation accumulates errors on the "high indexes" components resulting in a poor accuracy in terms of pixels values. TraIn and IDW instead interpolates globally the PSFs in the sense that the interpolation is not component-wise ; this implies less dependency of pixel-wise accuracy on the number of control points and in particular a better accuracy compared to the RBF interpolation where a few control points are used. However, this does not hold for the ellipticity. To understand these changes of performances, one needs to consider the shape parameters sensitivity relatively to the principal components. Specifically, for each PSF X i treated as an N l × N c matrix, we calculate the derivative of each ellipticity parameter along the lines passing though X i and parallel to each of the principal components P j also treated as matrices :

de k (X i +tP j ) dt , k = 1, 2.
The analytic expressions can be found in Appendix H. For each principal component P j , we denote disp j the dispersion of the PSFs set projected on the line passing through 0 and directed by P j . We define the sensitivity of the ellipticity component e k , k = 1, 2 with respect to the principal component P j as follows :

V(e k , P j ) = disp j (D + K) i | de k (X i + tP j )(0) dt |. (VI.34)
It measures how much small perturbations along each of the principal component is susceptible to make the ellipticity vary. In other words, it measures the sensitivity of the ellipticity to errors on the principal components coefficients. The sensitivities for the 40 first principal components are plotted in Fig.

VI.9 for the two ellipticity components. It is interesting to see that the plots are not monotonically decreasing. We can draw an important remark from this observation : the most important features in terms of pixel-wise error are not necessarily the most influential in terms of shape. The ellipticity parameters have average magnitudes of order 10 -2 . Therefore, one can tell from Fig. VI.9 that the ellipticity is significantly sensitive to only a few principal components within the 10 first one. Therefore, as long as the coefficient-wise interpolation is accurate for these low indexes principal components while staying in a reasonable range in general, the final result is accurate in terms of ellipticity. This explains why the RBF interpolation maintains a good accuracy on the ellipticity.

As to the size, it is clearly determined by the brightest structures on the PSF, namely the main lobe and the first brightest ring. The Figure VI.10 shows that the brightest ring information is distributed onto several components, in particular the first and the second one. One expects a bias in the RBF interpolation from processing this information components-wise. On the contrary, the brightest features are globally and more accurately modeled through the barycentric coordinates calculation step, since they influence it the most. Therefore, TraIn gives the most accurate results in terms of size. manifold, this parameter can be set to the manifold intrinsic dimension which is 2 in this example ; however this is tricky since the ideal density depends on the manifold geometric complexity which is a priori not known. We define the "extrinsic dimensionality" d ext of the data set as the dimension of the smaller subspace containing the data set. This can be approximately estimated using a PCA. Then we set d to min(p, d ext ).

F Practical considerations

Number of neighbors

The proposed method gives quite stable results with respect to the number of neighbors. However, Fig.

VI.11 suggests a potential improvement of the neighbors selection. For a given position of interpolation, the selection of neighbors can be pictured as follows :

-one initializes a null radius sphere centered on the interpolation position in the fov ;

-one increases the sphere radius until it includes a given number of observations locations.

Yet, Fig. VI.11 shows that the ellipticity parameters changes faster in certain directions than in others. This suggests to grow the neighborhoods according to the ellipticity parameters spatial gradients in such a way to include more neighbors from directions where the ellipticity parameters vary slowly.

Ground metric weight parameter We consider the parameter β in Eq. VI.13. β is the cost of matching two neighbors pixels (assuming a 4connected neighborhood) with identical intensities. Let consider two images x = (x i ) 1≤i≤N and y = (y i ) 1≤i≤N and the associated point clouds X and Y in the sense of Section C.

3. If ∀i ∈ 1, N |x i -y i | ≤ β then W 2 (µ X , µ Y ) =
xy 2 . Thus, β somehow determines how much the Wasserstein geometry is susceptible to depart from the pixels domain euclidean geometry. For sufficiently small neighborhoods around a given location in the fov, the pixels domain l 2 norm describes accurately the local geometry of the PSFs manifold. Therefore, β has to be chosen so that the Wassertein metric is equivalent to the pixels domain l 2 norm in any small vicinity in the fov. Hence the following procedure for setting γ :

-determine the two closest PSFs in the fov x i and x j ; -set γ to x i -x j ∞ .

F.2 Transportation issues

We recall the definition of the sliced Wasserstein distance between two clouds X and Y :

SW 2 (µ X , µ Y ) 2 = S d-1 W 2 (µ Xu , µ Yu ) 2 du, (VI.35) where S d-1 = {x ∈ R d / x = 1}, X u = {u T x i , i = 1...N } ⊂ R N and Y u is similarly defined.
Pairwise distances computation : assignments discrepancies The continuous integration being intractable, in practice, the sliced Wasserstein distance is calculated as a discrete sum :

SW 2 (µ X , µ Y ) 2 = u i ∈Ω⊂S d-1 W 2 (µ Xu i , µ Yu i ) 2 .
(VI.36)

For each vector u i we recall that the quadratic Wasserstein distance takes the form

W 2 (µ Xu i , µ Yu i ) 2 = N j=1 ((u T i X)[j] -(u T i Y)[σ * u i (j)]) 2 , (VI.37)
where σ u i is a permutation of 1, N . We use the stochastic gradient descent algorithm proposed in [Rabin et al. 2011] for calculating SW 2 . The algorithm estimates a local minimum Y * of the functional J Y (Z) = SW 2 (Z, Y) 2 in the vicinity of X and the approximated Wasserstein distance is computed as

W 2 (X, Y) ≈ X -Y * 2 . (VI.38)
At each gradient step, the "sliced assignment" (σ * u i ) i are updated. To guarantee the convergence, one can use a decreasing step size ≡ 1 n a for a ∈ 1/2, 1 , n being the iteration index [START_REF] Bottou | Online learning and stochastic approximations[END_REF]].

The algorithm succeeds in computing an assignment between the clouds X and Y if the slices assignments are identical at convergence and if the stationary point Y * represents the same point cloud as Y ; these two conditions are summarized below :

SW 2 (µ Y * , µ Y ) 2 = u i ∈Ω⊂S d-1 N j=1 ((u T i Y * )[j] -(u T i Y)[σ * u i (j)]) 2 = 0 and (VI.39) for u i = u j , σ * u i = σ * u j . (VI.40)
The displacement interpolation between the clouds represented by X and Y is realized by performing a linear interpolation between the matrices X and Y * . We define the discrepancy support as the following set

D = {k ∈ 1, N /∃(u i , u j ) ∈ Ω 2 /σ * u i [k] = σ * u j [k]} (VI.41)
and define the assignment discrepancy as as the support size |D|.

As observed in [Rabin et al. 2011] and [Bonneel et al. 2015], we do not found it necessary in practice to use a decreasing step size. We observe as well that the stationary points seem to always satisfy |D| = 0. However, it is not rare that the sequence of iterates oscillates around a point which has a nonzero discrepancy support size in which case, the algorithm takes a potentially long time to reach a stationary point satisfying |D| = 0.

In case |D| > 0 for the final iterate Y * , it mixes information from different pixels in the image associated to Y. This translates into visual artifacts when Y * is transformed into an image (see VI.15). One can see such an example in Fig. VI.12. This is a troublesome point since in the considered astronomical application, the systematic errors due to the PSFs estimation constitutes one of the bottlenecks.

Besides, the algorithm is stochastic and energy minimized is non-convex ; therefore for the same input point clouds, the algorithm converges to a different point at each run. We propose to reduce the final assignment discrepancy with a fixed number of iterations and the stationary point variability by improving the gradient descent initialization. Taking the two point clouds X and Y, we extract from each cloud the points corresponding to the pixels comprised in a small rectangular window (typically of size 20 × 20) around the corresponding images centroids (see Fig. VI.13). Let denote ind X and ind Y the sets of corresponding indexes and ind X the complement of ind X in 1, N . Then we find the optimal assignment between X[ind X ] and Y[ind Y ] using the Hungarian algorithm [Kuhn 1955], which is that fastest known procedure for solving exactly the assignment problem. However, it has a roughly cubic complexity in its improved version, which restricts its practical use to small scale problems. We denote σ the optimal mapping : Then we initialize the stochastic gradient descent method aforementioned with the matrix X 0 defined as follows :

X[ind X [i]] → Y[ind Y [σ(i)]], ∀i ∈ 1, |ind X | . (VI.42)
X 0 [:, i] = X[:, i] if i ∈ ind X and X 0 [:, ind X [i]] = Y[ind Y [σ(i)]] for i ∈ 1, |ind X | .
(VI.43)

To quantify the impact of this pre-assignment, we compare the average time of execution of the algorithm, for 3 different pair of images, for 100 runs in each case, with and without optimizing the initialization. These times account for the hungarian algorithm when used. We set a maximum number of iterations of 20000. The discrepancies support size |D| might be non-zero when this number of iterations is reached. Therefore, we also compute the average of the final discrepancies support sizes. The number of projection directions was set to 30 and the initialization window size to 20 × 20. As for all the experiments presented, we used 42 × 42 PSFs images, which gives clouds of 1724 points each one to compare. The result is displayed in Fig.

VI.14.

This results follows from the simple observation that the higher is |D| for a given iterate, the more the stochastic gradient descent is oscillatory. The Figure VI.14 shows that it is possible to set the exact matching window size so that the Hungarian algorithm brings more benefits in terms of discrepancies reduction than computational complexity.

It is worth noting that the assignment obtained from the sliced Wasserstein transport depends on the number of directions and slices. As outlined Figure VI.14 -Sliced Wasserstein transport algorithm average performances ; x axis : average execution times in seconds ; y axis : average final discrepancies support size ; the red, green and yellow corresponds to the pairs of images from the top to the bottom respectively ; the triangle-shaped points corresponds the gradient descents with an optimized initialization ; the points are obtained by averaging these performances criteria over 100 runs for the 3 pairs of 42 × 42 images ; the proposed initialization dramatically reduces the time needed to compute an assignment between our images related point clouds with an average execution time ranging from 10 to 40s rather than 100 to 400s.

in [Bonneel et al. 2015], the more directions, the smoother is the assignment. Unsurprisingly, this dependency is reduced when involving the Hungarian algorithm. Actually in that case, we obtain exactly the same assignment independently of the number of directions, for two given PSFs images and an initialization window of size 20 × 20.

Velocity constrained displacement interpolation The Wasserstein barycenter of two point clouds is simply calculated by realizing a displacement interpolation (see Section C.2). Considering 2 PSFs, a simplified interpolation process in this framework would go as follows :

- In these example, the brightest PSF's ring is more elongated horizontally for the initial PSF and vertically for the final PSF. But as the displacement interpolation constraints particles to move along straight lines, the interpolated point clouds are shrunken compared to the extreme ones and consequently, the interpolated images have a more concentrated energy. Quantitatively, the interpolated PSFs tend to have l 2 /l 1 norms ratio that significantly exceed the range observed over the data set -note that the PSFs l 1 norms are constant and equal to 1. This problem can be approached from the ground metric angle ; precisely, one can think of choosing a ground metric that admits curved geodesics, making a more "compact" particles advection possible. However, the complex variation of the PSFs across the fov makes the choice of a better fitted ground metric non trivial. The displacement interpolation generates a sequence of matrices

(X i ) 1≤i≤T so that X 1 = X, X T = Y * and for i ∈ 1, T -1 , X i+1 -X i = µ i (Y * -X) for a step µ i .
In order to better preserve the point clouds densities in the displacement interpolation, we impose the points to have parallel velocity vectors at each time. To do so, we modify the displacement interpolation as follows :

1. INPUT : a maximum step size µ max 2. Initialization : X 0 = X, i = 0.

. WHILE i ≤ i max : V i = Y * -X i Compute an eigenvector u i corresponding to the highest eigenvalue of V i V T i V i = u i u T i V i / u i 2 2 Line search : µ opt = argmin µ X i + µ V i -Y * 2 2 = -V i ,X i -Y * V i 2 2 X i+1 = X i + min(µ opt , µ max ) V i i = i + 1.
In the above procedure, all the particles are moved in parallel with u i at each iteration. Moreover, X i+1 -Y * 2 2 ≤ X i -Y * 2 2 . Therefore the procedure moves in fact the initial point cloud toward the final one. We do not discuss how close the sequence generated can come to Y * . However, in the tests we performed, the distance X i -Y * always gets down to the machine numerical precision. In Fig. VI.15, the time parameter is defined as

t(i) = 1 - X i -Y * 2 X -Y * 2 , (VI.44)
which makes sense in a dynamic model of constant speed advection along straight lines. However, in the modified displacement interpolation, the matrices (X i ) 1≤i≤imax does not necessarily follow a linear path in the matrix space. To make a direct comparison possible between the two displacement interpolations, the time parameter definition has to be generalized in order to account for a possible curvature of the path. We do so by involving the length of the curve obtained by joining each matrix of the sequence to the following one with a straight line :

t(i) = 1 - imax j=i X i -X i+1 2 imax j=0 X i -X i+1 2 , (VI.45)
with the convention X imax+1 = Y * . This definition generalizes VI.44. We compare the two displacement interpolations in Fig. VI.16. The parallel velocity constraint quantitatively reduces the shrinkage effect since the interpolated images l 2 /l 1 norms ratios does not exceed those of the reference images as much as with the regular displacement interpolation. Not surprisingly, this constraint yields smoother barycenters. This constraint can be integrated directly into the sliced OT algorithm. However, we observe that this considerably increases the number of iterations needed to converge and most importantly, the generated sequence mostly converges to a solution which is not a global minimum of the functional J Y (Z) = SW 2 (Z, Y) 2 . We note that this constraint can be naturally integrated into the framework adopted in [Papadakis et al. 2013] since velocities are directly manipulated. However, the eulerian discretization would require manipulating 4 dimensional arrays, making the approach intractable. This constraint has not been used for generating the plots shown in Section E. Indeed, for that experiment, we considered a density of known PSFs that makes the PSFs variations in a given neighborhood in the fov very smooth in which case the shrinkage effect disappears.

Wasserstein barycenters computation

As explained in Section D, the interpolated PSF's computation involves calculating a sliced Wasserstein barycenter. Let consider a set of point clouds (Y i ) 1≤i≤p . As for the pairwise Wasserstein distances, an approximation of their Wasserstein barycenter is approximately calculated by finding a local stationary point of the functional Rabin et al. 2011], for some barycentric weights (w i ) 1≤i≤p . The higher is p, the more degenerated is this functional because of its combinatorial nature, which increases the assignments discrepancies. But as mentioned in Section C.2.2, there is no efficient method for calculating exactly a Wasserstein barycenter of more than two point clouds, even for small scale problems. Therefore, the previously described strategy for speeding the convergence up can not be directly extended to this case. We propose to approximate the Wasserstein barycenter of a set of more than two point clouds by computing a sequence of "2 points" Wasserstein barycenters. This relies on the local isometry assumption made in Section D.1. Indeed, the barycenter of a set of vectors (r i ) 1≤i≤p associated with the barycentric weights (w i ) 1≤i≤p can be calculated through the following procedure : norms ratios are flattened compared to those obtained with the regular displacement interpolation ; moreover, the interpolated images are smoother.

J Y 1 ,...,Yp (Z) = p i=1 w i SW 2 (Z, Y i ) 2 [
w = w + w i i = i + 1.
Indeed, one can check that at the i th iteration r bar is updated to i k=1 w k i j=1 w j r k . If the Wasserstein space is isometric to a euclidean space in the neighborhood of the point clouds involved in the barycenter's calculation, then one can apply a similar procedure in the Wasserstein space to calculate the barycenter, hence the following scheme :

1. Initialization :

Z bar = Z 1 ; i = 2, w = w 1 2. WHILE i ≤ p : Z bar = argmin Z w i SW 2 (Z, Y i ) 2 + wSW 2 (Z, Z bar ) 2 w = w + w i i = i + 1.
This way we can take advantage of the previously accurate initializing proposed in the Wasserstein barycenter approximation and make use of the 1D displacement interpolation. The computed barycenter slightly changes depending on the ordering of the clouds Y i . We choose the ordering so that

w 1 ≥ • • • ≥ w p .

G Reproducible research

In the spirit of participating in reproducible research, the data and the codes used to generate the plots presented in this paper will be made available at http://www.cosmostat.org/software/.

H Conclusion

We introduce TraIn (Transport Interpolation) which is a data field interpolation method based on Optimal Transport and making use of some Manifold Learning ideas. We consider the interpolation of a PSF field. The local geometry of the PSF field is characterized using approximated Wasserstein distances. From these, we derive low dimensional local euclidean embeddings of the PSF field which is then mapped to the instrument field-of-view using a thin-plate interpolation. This mapping gives one the embedded coordinates of the unknown PSFs, from which a representation of these PSFs in the quadratic Wasserstein metric is determined. Finally, the interpolated PSFs are calculated through nested displacement interpolations.

We compared TraIn to the Inverse Distance Weighting method and to a component-wise thin-plate interpolation of PCA coefficients. The tests were made on a set of realistic Euclid-like PSFs. We show that the proposed method is globally the most accurate in terms of pixel domain error and shape. In particular, TraIn is in average several order of magnitudes more accurate than the two other methods in terms of pixels errors when a few number of "neighbors PSFs" are used for the interpolation ; it is therefore more robust to the PSF field sampling. We also introduce a velocity constraint displacement interpolation for mitigating the unnatural shape shrinkage that might occur when interpolating objects with different major axis using an euclidean ground metric.

A natural extension of this study would be to compare these interpolation methods in a setting where the reference PSFs are not perfectly known.

Chapitre VII

Conclusion

A Methods and results

This thesis addressed the problem of estimating the Point Spread Function across an optical instrument field-of-view (fov) solely from unresolved images of compact objects such as stars.

In Chapter IV, we proposed a method for estimating a well sampled PSF from aliased, randomly shifted and noisy measurements of the same PSF. We show that using sparsity and positivity constraints yield a significantly better preservation of the PSF shape over popular methods in the astronomy community, especially at low SNRs.

In Chapter V, we move toward a more realistic setting. Indeed, we consider the problem of jointly estimating a set of different PSFs from only one measurement for each PSF. Thus, we introduce a dimension reduction and super-resolution method which, additionally to the priors used in the multiple frame setting, takes advantage of the PSFs field compressibility. This, again, yields significant improvements in terms of the PSFs shapes at low SNR, with respect to existing methods.

In Chapter VI, we consider the problem of estimating the PSF at an arbitrary location in the instrument fov, given the PSFs on a finite set of locations. This amounts to estimating the mapping between the fov coordinates space and the PSFs underlying manifold. To that end, we propose an interpolation framework which firstly derives local low dimensional representations of the training data based on approximated pairwise quadratic Wasserstein distances. Assuming that this unfolds the PSFs manifold locally, the low dimensional representations of the unknown PSFs can be estimated using standard interpolation methods. From these, the barycentric coordinates of the unknown PSFs in the Wasserstein metric are derived and the interpolated PSFs are calculated as (approximated) Wasserstein barycenters. This approach yielded a remarkable accuracy in terms of the PSFs shapes and pixels errors. However, it is computationally demanding.

We emphasize that the methods developed are non parametric and can be used for other optical instruments, although the numerical experiments were focused on the Euclid telescope PSFs. The associated codes will be made available at http://www.cosmostat.org/software/.

B Perspectives

Several perspectives follow naturally from this work and applications in astronomy and beyond.

In astronomy : The interpolation problem can be alternately addressed as an extension of RCA (see Chapter V) by solving a larger inverse problem involving a "semi-blind" deconvolution of the galaxies present in the same fov. We recall that RCA computes a matrix S which columns span a subspace embedding the PSFs manifold and a matrix A of representation coefficients of the PSFs at the positions of the unresolved stars images present in the fov. The ultimate aim is to estimate the PSFs at the galaxies locations. Assuming that these PSFs also belong the S's columns space, one only need to determine their representations coefficients, which we denote by a matrix A g . Given an estimate of A, A g is strongly constrained, due to the PSFs field's regularity and the PSFs positivity prior. Given an estimate of A g and S, one can restore the galaxies. Finally, given an estimate of the galaxies, one can refine A g . Thus, it is easy to derive a global scheme which both restores and interpolates the PSFs and deconvolves the galaxies.

More generally : The methods proposed have to be extended to account for the wavelength dependency of the PSFs. Specifically, an unresolved image of a compact light source can be modeled as y = n i=1 F(a i x i ) + n, (VII.1)

where each x i is the PSF at the position of object in a narrow wavelength band centered on a wavelength λ i . The vector a = (a i ) 1≤i≤n is the discretized emission spectrum of the light source. n has to be sufficiently large for the PSF to be approximately constant on each band. n is the noise and F is a linear degradation operator. Let assume that F is the identity. Calculating the monochromatic PSFs x i from y is an ill-posed inverse problem, even knowing a. However, this problem is highly structured. Indeed, assuming that λ 1 < • • • < λ n , then for i > 1, the PSF x i is approximately obtained by rescaling the PSF x i-1 as illustrated in Fig. VII.1. This follows from the dilatation property of the Fourier transform ; we refer to the introductory chapter for the analytical form of a monochromatic PSF. It is not exactly a rescaling because of the wavelength dependency of the wavefront error ; again, we refer the reader to Chapter I. However, we can reasonably assume that in a continuous setting, the PSF x i-1 can be transformed into the PSF x i through a Figure VII.1 -Simulated Euclid optical PSFs at increasing wavelengths from the left to the right between 500 and 900 nm monotonic rearrangement. Thus, rather than computing directly the x i s, one can approach the problem by computing a sparse and positive transportation plan between x 1 and x n , which generates the intermediary PSFs by displacement interpolation and approximately performs a monotonic rearrangement of x 1 into x n . One can think of different tools to enforce the monotonicity, among which the helmoltz-hodge decomposition, in the same spirit as [START_REF] Zhu | [END_REF]]. Conveniently, this approach would implicitly integrates the mass conservation constraint between the PSFs. Furthermore, as seen in the introductory chapter, more physical priors could be involved in the PSF's restoration. A simple example is the fact that a PSF's Fourier transform is an even function, which can be translated into a linear constraint. A bolder approach would be to turn the PSF restoration into a phase recovery problem, and use prior knowledge on the phase-shifting function, like for instance its relationship with Zernike polynomials. This is particularly sound if out-of-focus images are available, which can be easily obtained in microscopy applications [START_REF] Bostan | [END_REF]].

Teaching

-Scientific activity tutor in primary schools for the Center of Research and Inter-disciplinarity (CRI), 2013-2015.

-Teaching assistant at ENSAE ParisTech in Numerical Analysis, 2015-2016.

Miscellaneous

-2nd price for the "Concours découvrir l'Univers" of the SF2A with the CRI, May 2014.

Annexe B In this appendix, we give the general convex analysis material relevant to our work. We consider a finite-dimensional Hilbert space H equipped with the inner product ., . and associated with the norm . . Let H be a finitedimensional Hilbert space (typically a real vector space) equipped with the inner product ., . and associated with the norm . . A real-valued function F defined on H is -proper if its domain, as defined by dom F = {x ∈ H/F(x) < +∞}, is non-empty ;

Convex analysis

-lower semicontinuous (LSC) if lim inf x→x 0 F(x) ≥ F(x 0 ).

We define Γ 0 (H) as the class of all proper LSC convex real-valued function defined on H.

A Proximal calculus

Moreau (1962) introduced the notion of proximity operator as a generalization of a convex projection operator. Let F ∈ Γ 0 (H). Then the function y → 1 2 αy 2 + F(y) achieves its minimum at a unique point denoted by prox F (α), (∀α ∈ H). The operator prox F is the proximity operator of F. It is clear from the definitions that the proximity operator of ι C is the orthogonal projector onto C. Thus, for C = S x (0) , which is defined in and F 2 would constrain this solution based on prior knowledges. It has been shown [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] that if F 1 is differentiable with a β-Lipschitz continuous gradient, then the problem (B.6) admits at least one solution and that its solutions, for γ > 0, verify the fixed point equation,

x = prox γF 2 (x -γ∇F 1 (x)).

(B.7)

This suggests the following iterative scheme,

x n+1 = prox γnF 2 (x n -γ n ∇F 1 (x n )), (B.8)
for appropriate values of the parameter γ n . This type of scheme is known as forward-backward (FB) algorithm : a forward gradient step using F 1 and a backward step involving F 2 through its proximity operator. Some variants of FB algorithms that have been shown to converge to a solution of (B.6) can be found in [START_REF] Bauschke | [END_REF].

One may refer to [START_REF] Bauschke | [END_REF] for other proximity operator properties and examples.

B Convex conjugate

Definition : let F be a real-valued function defined on H. The function F * : y → max x x, y -F(x) is the convex conjugate of F ; it is also known as the Legendre-Fenchel transformation of F.

Properties :

-Moreau identity : for F ∈ Γ 0 (H) and λ ∈ R * + , prox λF (x) + λprox 1 λ F * ( x λ ) = x ; -Fenchel -Moreau theorem : if F ∈ Γ 0 (H), F = F * * .

Annexe C

Multiple frame super-resolution

A First guess noise

We keep the same notations as in Section B.3.2 in Chapter IV. We note x wi as the i th scale of x (0) n wavelet transform in W. As the noise in x (0) n is correlated, the estimation of σ i is not straightforward, so we proceed as follows :

1. σ 0 i = 1.4826 MAD(x wi ) 2. xwi = argmin The factor 1.4826 comes from the assumption that the noise is approximately Gaussian. Finally, the factor k must be sufficiently high so that all the noise will remain in the residual x wi -x wi . We took k = 5. The problem in step 2 has a closed-form solution. For the finer scales, σ 0 i is quite close to σ i . But for coarser scales, σ 0 i is significantly overestimated (it might be more than 10 times greater than σ i ). We implicitly assumed that the noise in the wavelet scales is stationary, which is reasonable apart from the edges effects due to the wavelet transform.

B Positivity

We can drop the positivity constraint by simply solving

C Reweighting

As stated in Chapter IV, Section B.1, the reweighting scheme used in SPRITE is meant to mitigate the bias due to the l 1 norm penalty. To verify this, we basically did the same as for the positivity in the previous section. Thus, we ran Algorithm IV.2 with K max = 1 and K max = 2 and we compute in each case the mean correlation coefficient in Pearson sense [Rodgers & Nicewander 1988] This appendix details the practical resolution of optimization problems involved in the RCA algorithm. F for some bounded linear operator M. F ∈ Γ 0 (R n×r ), G 1 ∈ Γ 0 (R n×p ) and L 1 is a bounded linear operator. Moreover, H is convex, differentiable and has a continuous and Lipschitz gradient. This problem can be solved efficiently using the primal dual algorithms introduced in [Combettes et al. 2014] for instance. One only need to be able to compute λF and αG * 1 proximity operators, for some given positive reals λ and α and H's gradient :

A Components estimation problem

-prox λF (S) = (ŝ ij )1≤i≤n For an arbitrary dictionary Φ s , we instead consider the following generic formulation of the problem :

B Coefficients estimation

We consider the step 8 in the Alg. V.1 in Chapter V. The problem takes the generic form : where J is convex, differentiable and has a continuous and Lipschitz gradient and α ∈ R r×q . This problem is combinatorial and its feasible set is nonconvex. For typical data sizes in image processing applications and tractable processing time, one can at best reach a "good" local optimum. There is an extensive literature on optimization problems involving the l 0 pseudo-norm.

We propose an heuristic based on quite common ideas now and which appears to be convenient from a practical point of view. Let α * be a global minimum of Problem D.3. For a vector M ∈ R r×q , we define its support as Supp(M) = {(i, j) ∈ 1, r × 1, q /|M[i, j]| ≥ 0}. (D.4)

We note E α * the set of r × q real matrices sharing the support of α * : E α * = {M ∈ R r×q /Supp(M) = Supp(α * )}.

(D.5) E α * is a vector space. In particular, E α * is a convex set. Therefore, α * is a solution of the following problem : for a given subspace E, which is a convex problem. This scheme can be viewed as an iterative hard thresholding [Blumensath & Davies 2008], with a decreasing threshold [START_REF] Mancera | [END_REF]]. Yet, it is quite easy to get an upper bound of the support size -related to the parameters η l in Problem D.3 -from the data. Depending on the time one is willing to spend on the coefficients computation, this yields convenient choices for the function f .

U k = β k -ρ -1 ∇J (β k ) 5: α k+1 = proj E f (k) ,U k (U k ) 6: t k+1 = 1+ √ 4t 2 k +1 2 7: λ k = 1 + t k -1 t k+1 8: β k+1 = α k + λ k (α k+1 -α k ) 9:

Notch filter approximation

In this appendix, we explain why the functional Ψ e,a introduced in Chapter V, subsection C.2 can be approximated with the functional Ψ e,a . We reuse the same notations. We consider the 1D case. The samples (u i ) 1≤i≤p are uniformly spaced scalar. We assume that u 1 < • • • < u p . We note ∆ = u 2 -u 1 . Thus, Given the symmetry of ψ e,a with respect to k + 1, we focus on the above difference for j ≤ k. We further assume that j = 1. Then, Eq.E.6 simplifies to We see that the higher is k (we recall that p = 2 * k + 1) and the closer j is to k, the smaller is the error. Therefore, we use t v as an approximation for h, up to boundaries errors.

ψ i = ψ p-i+1 = -1 |k + 1 -i| e

Annexe H

Ellipticity parameters directional derivatives

The ellipticity parameters used in the Chapters IV to VI can be rewritten in the following way : . We derive the following expressions : where c t =< U 4 , X i + tP j >< U 3 , X i + tP j > -< U 1 , X i + tP j > 2 -< U 2 , X i + tP j > 2 a 1 =< U 5 , X i >< U 3 , P j > + < U 3 , X i >< U 5 , P j > -2(< U 1 , X i >< U 1 , P j > -< U 2 , X i >< U 2 , P j >)

e 1 (X i ) = < X i , U 5 >< X i , U 3 > -< X i , U 1 > 2 + < X i , U 2 > 2 < X i , U 4 >< X i , U 3 > -< X i , U 1 > 2 -< X i , U 2 > 2 (H.1) e 2 (X i ) = 2(< X i , U 6 >< X i , U 3 > -< X i , U 1 >< X i , U 2 >) < X i , U 4 >< X i , U 3 > -< X i , U 1 > 2 -< X i , U 2 > 2 , ( H
de 1 (X i + tP j ) dt =
a 2 = 2(< U 5 , P j >< U 3 , P j > -< U 1 , P j > 2 + < U 2 , P j > 2 )

b 1 = 2(< U 6 , X i >< U 3 , P j > + < U 3 , X i >< U 6 , P j > -< U 1 , X i >< U 2 , P j > -< U 2 , X i >< U 1 , P j >)

b 2 = 4(< U 6 , P j >< U 3 , P j > -< U 2 , P j >< U 1 , P j >)

d 1 =< U 4 , X i >< U 3 , P j > + < U 3 , X i >< U 4 , P j > -2(< U 1 , X i >< U 1 , P j > + < U 2 , X i >< U 2 , P j >)

d 2 = 2(< U 4 , P j >< U 3 , P j > -< U 1 , P j > 2 -< U 2 , P j > 2 ).

-

  CG : Conjugate Gradient -CMB : Cosmic Microwave Background -CTI : Charge Transfert Inefficiency -ESA : European spatial agency fov : field-of-view -FWHM : Full Width at Half Maximum i.e. : id est -GMRA : Geometric Multi-Resolution Analysis -HLLE : Hessian-based Local Linear Embedding --w.r.t. : with respect to

  I.1).

Figure

  Figure I.1 -The Universe content.

  Fig. .

Figure I. 2 -

 2 Figure I.2 -Weak gravitational lensing.

Figure I. 3 -

 3 Figure I.3 -Sun-Earth system Lagrange points. Credit : nasa.gov

Figure I. 4 -

 4 Figure I.4 -Planar projection of the sky : the distortion associated to thisprojection with respect to the real world 3D objects is not considered.

  I.5) is close to the wavelength.

Figure I. 5 -

 5 Figure I.5 -Diffraction : Huygens-Fresnel principle.

Figure I. 6 -

 6 Figure I.6 -Diffraction : Huygens-Fresnel principle.

  Figure I.7 -"Black box" linear optical system ; U i : input wave ; U o : output wave

Figure

  Figure I.9 -"Appeareance manifold" : all the images are generated with a single facial pose, moving a single light source.

Figure

  Figure I.10 -S-curve data set.

A. 1

 1 Figure II.1 -For instance natural images are in general piece-wise smooth ; the bottom panel shows an highly structured speech signal which can be represented using a few parameters which in fact characterize the vocal organ's state during the phonation.

  is the degradation operator, b represents the noise and/or the modeling error, -• is some composition operator in R p .

Figure

  Figure II.2 -Deconvolution by direct convolution kernel inversion ; from the left to the right : reference undistorted image and blurring kernel, blurred image, restored images by direct inversion in Fourier domain, for snr = 10 12 , 10 9 and 10 6 respectively.

Figure

  Figure II.3 -Ill-posedness and ill-conditioning in inverse problems : different signals in the undistorted signals space can be mapped to the same observation, which implies ill-posedness ; more generally, significantly different signals in the undistorted signals space can be mapped to very close observations as a consequence of ill-conditioning.

1

  cost : C(y, Mx) = y -Mx 1 [Yang et al. 2009]. -Similarly, if the signal is dominated with the so-called Poissonian noise, the natural cost function in a bayesian sense is the Kullback-Leibler divergence : C(y, Mx) = p i=1 y[i] log( [

  26) for i ∈ 1, t + . The weights variations as a function of the singular values is illustrated in Fig. II.4, for different numbers of iterations. For k > 1, these functions have a step-like shape with an inflection point localized at

Figure II. 4 -

 4 Figure II.4 -Curve of the weight in Eq.II.26 as a function of the singular values, for different number of iterations.

Figure

  Figure II.5 -On the left, the first detected gravitational waves signal, measured by the LIGO Hanford (credit : www.ligo.org) ; the x axis is the time and the y axis is the gravitational waves strain, which quantifies the space-time distortion. On the right, the DCT of this signal (truncated to the first 500 hundred samples for clarity, the discarded values being negligible).

Figure II. 6 -

 6 Figure II.6 -The sorted magnitude of Fig. II.5's signals. The gravitational is clearly sparser in the sense of Eq.II.31 in the DCT domain, although more targeted transforms might have been used (see for example the chirplet transform in [Mann & Haykin 1995]).

  Figure III.1 -S-curve data set : Isomap estimates the riemanian distances based on a nearest neighbor graph ; when the samples are sufficiently dense, this prevents one from bypassing the manifold's geometry.

Figure

  Figure III.2 -Minimal effort packaging.

  III.21) and we note ν = T #µ (III.22) T is then referred as a transport map between µ and ν. This notion is illustrated in Fig.III.3.

Figure

  Figure III.3 -Measure pushforward : the shaded surfaces should have the same area.

Y

  dΠ(x, y) = dµ(x) and X dΠ(x, y) = dν(x). (III.24) dΠ(x, y) represents the elementary mass transferred from x to y. From this definition, the Kantorovich Optimal Transport problem can be stated as follows : min T ∈P(µ,ν) X×Y c(x, y)dΠ(x, y), (III.25) where P(µ, ν) is the set of transportation plans between µ and ν. Note that the probability distribution of X × Y µ ⊗ ν : (x, y) → µ(x)ν(y) (III.26)

Figure

  Figure III.5 -Three gaussian probability measures and their euclidian and quadratic Wasserstein barycenters.

  Fig. IV.1.

Figure

  Figure IV.1 -General observation model. See Section A.2 for a detailed description

Figure

  Figure IV.2 -Adapted observation model. This time we consider an isolated star as an input, and unlike in the general model, the output of the block "Blur" is the PSF.

Figure

  Figure IV.3 -SR general scheme. The deblurring does not apply to our case, since we want to precisely estimate the blur.

Figure IV. 4 -

 4 Figure IV.4 -Simulated optical point-spread function.

  IV.5 and Fig. IV.6 below). These LR images are of size 84×84. Different levels of WGN are added.

Figure

  Figure IV.5 -Critically sampled PSF.

Figure

  Figure IV.6 -PSF sampled at Euclid resolution with different offsets and noise.

  (a) Super-resolved PSFs. (b) Absolute Error Maps.

Figure

  Figure IV.7 -PSF reconstruction and error map at 30dB for three methods : from the left to the right, shift-and-add, PSFEx, and SPRITE. The error image standard deviation is at least 30% smaller with SPRITE.

( a )

 a Mean absolute error on the first ellipticity parameter. (b) Standard deviation of the absolute error on the first ellipticity parameter. (c) Mean absolute error on the second ellipticity parameter. (d) Standard deviation of the absolute error on the second ellipticity parameter.

Figure

  Figure IV.8 -Errors in log on ellipticity parameters versus the SNR. For SNR=10dB, SPRITE achieves around 6dB less than others methods, which corresponds to a factor of e 6 on a linear scale.

Figure

  Figure IV.9 -Mean absolute error on the full width at half maximum (FWHM) in percent. SPRITE achieves on average an error of 6% on the FWHM which is 2% less than PSFEX in average.

  r, (V.9) where (υ l ) 1≤l≤N is a set of vectors spanning R p . Equivalently, we can write A = αV T , where V = [υ 1 , . . . , υ N ] and α is a r × N matrix (see Fig. V.1).

  (a) Direct domain samples. (b) Discrete Fourier Transform (DFT) entry-wise moduli.

Figure V. 2 -

 2 Figure V.2 -Notch filters examples for different values of the parameter e in Eq. V.10 and V.11. The parameter a is set to 1.

Figure

  Figure V.3 -Simulated PSFs distribution across the FOV.

( a )

 a Average error on the ellipticity vector.(b) Dispersion of the ellipticity vector.

Figure

  Figure V.4x axis : SNR (see section D.2) ; y axis : log 10 (E γ ) for the left plot, log 10 (B γ ) for the right plot.

( a )

 a Average absolute error on the size. (b) Dispersion of the errors on the size.

Figure

  Figure V.5x axis : SNR ; y axis : E S for the left plot, σ S for the right plot.

  V.3 at an SNR of 40. This assess the behavior of the algorithms with respect to spatial clustering and sparse data distribution. One can see in Fig. V.6 examples of simulated observed PSFs from different areas in the FOV. For each of these observed PSFs, the reconstructed PSFs for each method are shown in Fig. V.7.

Figure V. 6 -

 6 Figure V.6 -Input PSFs at different locations in the FOV for a SNR = 40. The corresponding reconstructed PSFs can be seen in Fig. V.7

Figure V. 7 -

 7 Figure V.7 -PSFs reconstructions : from the left to the right : original, GMRA, PCA, PSFEx, RCA ; from the bottom to the top : 2 "local" PSFs reconstructions, 2 corner PSFs reconstructions, 2 center PSFs reconstructions. The observed corresponding PSFs can be seen in Fig. V.6

( a )

 a Average error on the ellipticity vector.(b) Dispersion of the ellipticity vector.

Figure V. 8

 8 Figure V.8x axis : SNR (see section D.2) ; y axis : log 10 (E γ ) for the left plot, log 10 (B γ ) for the right plot.

( a )

 a Average absolute error on the size. (b) Dispersion of the errors on the size.

Figure

  Figure V.9x axis : SNR ; y axis : E S for the left plot, σ S for the right plot.

Figure V. 10 -

 10 Figure V.10 -Average normalized least square error

Figure V. 11 -

 11 Figure V.11 -Input PSFs at different locations in the field for a SNR = 40.

Figure

  Figure V.12 -PSFs reconstructions : from the left to the right : original, PSFEx, RCA ; from the bottom to the top : 2 "local" PSFs reconstructions, 2 corner PSFs reconstructions, 2 center PSFs reconstructions. The observed corresponding PSFs can be seen in Fig. V.11

Figure

  Figure VI.1 -Simulated Euclid telescope PSFs at different locations in the fov

Figure

  Figure VI.2 -PSF interpolation problem

  Algorithme VI.1 Transport Interpolation (TraIn)1: Inputs : K PSFs (x k ) 1≤k≤K , K observations locations (u k ) 1≤k≤K , D interpolation locations (v k ) 1≤k≤D ,number of neighbors p, local dimensionality d ≤ p 2: Compute the weighting parameter β (see Eq.VI.14) 3: Compute the p neighbors of each location v i in the set (u k ) 1≤k≤K ; this results in a collection of sets of indices (I(v i )) 1≤i≤D 4: Transform the PSFs into point clouds (see Eq. VI.12) : (x k ) 1≤k≤K (X k ) 1≤k≤K 5: Compute the approximated Wassertein 2 distances between pairs of point clouds (X l , X m ) /(l, m) ∈ I(v i ) 2 for some i ∈ 1, D 6: Pour i = 1 to D faire 7:

Figure

  Figure VI.3 -Simulated PSFs distribution across the fov ; the blue squares represents the 300 observed PSFs and the red circles the 250 PSFs to be interpolated ; the coordinates are in degrees.

  ) 2 /D ; -the average absolute error on the size : E S = D i=1 |S(Im i ) -S( Im i )|/D in pixels ; -the average normalized mean square error : NMSE =

Figure

  Figure VI.4 -Normalized mean squared error : Y-axis log 10 (NMSE), X-axis number of neighbors.

Figure VI. 5 -

 5 Figure VI.5 -Average error on the ellipticity : Y-axis E γ , X-axis number of neighbors.

Figure VI. 6 -

 6 Figure VI.6 -Average absolute error on the size : Y-axis E S , X-axis number of neighbors.

  |Im i -Im i | for some i ∈ 1, N inFig. VI.7, which correspond to interpolations with 15 neighbors. It shows indeed that TraIn typically yields a substantially lower residual.

Figure

  Figure VI.7 -Absolute error images for 15 neighbors ; from the left to the right, TraIn, RBF, IDW, original image ; from the top to the bottom, best, average and worst case in terms of the relative accuracy of TraIn with respect to the pixels mean square error.

Figure

  Figure VI.8 -Spatial distributions of the coefficients of the PSFs used for the experiments relative to 4 principal components calculated over the "learning" PSFs set ; the top panel corresponds to the first and second principal components and the coefficients scatter plot ; the bottom panel corresponds to the 6th and 7th principal components and the coefficients scatter plot as well.

F. 1

 1 ParametersLocal dimensionality We refer to the parameter d involved in the step 7 in Algorithm VI.1. For an ideally dense data set in the sense of the underlying

Figure

  Figure VI.9 -Ellipticity components sensitivity to the PCA atoms. Yaxis : log 10 (V(e k , P j ))(see Eq.VI.34) ; X axis : PCA atoms indexes j ; moderate errors on the PCA coefficients have a weak impact on the ellipticity, especially for the "high indexes" PCA components.

Figure VI. 10 -

 10 Figure VI.10 -First and second principal components ; the brightest ring is broken apart between the two components.

Figure

  Figure VI.11 -Spatial distributions of the first (left) and second ellipticity (right) parameters of the observed PSFs.

Figure VI. 12 -

 12 Figure VI.12 -Assignment failure : on the top plot, the sequence of images associated with each iterate are represented by the red dots and the image associated with the ideal solution is represented by the blue dot ; these images are shown in the plan spanned by the two first principal components calculated from the whole images set ; the iterates end up oscillating around a point in the middle circle ; the bottom left plot shows the ideal image and the bottom right plot shows the absolute difference between the final image and the ideal one ; one can see that there is some non-negligible artifacts.

Figure

  Figure VI.13 -Two PSFs examples in a logarithmic scale ; we first use the hungarian algorithm to map the two point clouds derived from the pixels located in the small red rectangles in each PSF respectively ; the result is used for initializing the sliced transport algorithm on the point clouds associated with the full images (see VI.12).

  the PSFs are first converted into two matrices X and Y representing two different point clouds ; -the point clouds are matched according to the previous paragraph ; this yields a matrix Y * representing the same point cloud as Y ; -the point clouds interpolation is realized by moving each point in the first point cloud toward the matched point in the second point cloud along a straight line and at a uniform speed on the cloud ; this is done by performing a linear interpolation of the matrices X and Y * ; -the interpolated PSFs are obtained by converting interpolated clouds at different "times" into images. An example is shown in Fig.VI.15.

Figure

  Figure VI.15 -Displacement interpolation : "shrinkage effect" ; the interpolated PSFs have a narrower energy distribution than the initial and final ones.

  3

2 Figure

 2 Figure VI.16 -Constrained displacement interpolation : the two PSFs panels show original PSFs on the top and the interpolated PSFs for t=0.5 are shown on the bottom left for the regular displacement interpolation and the bottom right for the velocity constrained displacement interpolation ; the right plots show the evolution of the l 2 /l 1 norms ratio of interpolated PSFs as functions of the time parameter in the two examples respectively ; the velocity constraint reduces the shrinkage effect since the corresponding l 2 /l 1 norms ratios are flattened compared to those obtained with the regular displacement interpolation ; moreover, the interpolated images are smoother.

  . . . . . . . . . . . . . . . . . . . . . . . 145 B Convex conjugate . . . . . . . . . . . . . . . . . . . . . . . . 147

-

  The indicator function of a closed convex subset C of H is the function defined on H by ι C (x) = 0, if x ∈ C +∞, otherwise. (B.1)

  . . . . . . . . . . . . . . . . . . . . . . . . . 149 B Positivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 C Reweighting . . . . . . . . . . . . . . . . . . . . . . . . . . . 151

  i = 1.4826 MAD( x wi -x wi ).

  ∆ + x (0) ) + κ w(k) λ Φ∆ 1 (C.1) at step 6 in Algorithm IV.2 in Chapter IV. We did a similar numerical experiment as the one presented in Section C in the same chapter to quantify the impact of this constraint. The comparison is given in Fig. C.1. One can see that the positivity constraint is actually important from the point view of the ellipticity parameters at low SNR. This result is illustrated in Fig. C.2. With the positivity constraint, the reconstruction is less influenced by the negative oscillations in the data, due to noise ; thus it yields a better robustness. Even if the negative residual values in the final PSF are 1000 order of magnitude smaller than the peak value, it is sufficient to considerably bias the ellipticity measurements.

( a )

 a Mean absolute error on the first ellipticity parameter. (b) Standard deviation of the absolute error on the first ellipticity parameter. (c) Mean absolute error on the second ellipticity parameter. (d) Standard deviation of the absolute error on the second ellipticity parameter.

Figure C. 1 -

 1 Figure C.1 -Errors in log scale on ellipticity parameters versus the SNR.The positivity constraint significantly improves the accuracy.

Figure C. 2 -

 2 Figure C.2 -Map of negative values in the PSF reconstruction map (in absolute value) at 15dB. On the left, SPRITE has a positivity constraint ; on the right, SPRITE without positivity constraint.

  between the reference images and the reconstructions for different SNR. The result is given in Fig. C.3. As expected, the reweighting improves the correlation and consequently, reduces the global bias on the reconstruction.

Figure C. 3 -

 3 Figure C.3 -Mean correlation coefficients (×100) between the SPRITE PSF reconstructions and the reference images versus the SNR ; the reweighting reduces the global bias

  We consider the step 5 in the Algorithm V.1 in Chapter V. If Φ s = I n , the problem of estimating components takes the following generic form :min S F(S) + G 1 (L 1 (S)) + H(S), (D.1) with F(S) = r i=1 w i s (c) i 1 , G 1 = ι R n×p + , L 1 (S) = SA and H(S) = 1 2 Y -M(S) 2

1≤j≤p,

  with ŝij = SoftThresh λw j [i] (s j [i]) ; -prox αG * 1 (Z) = Z -(Z) + -∇H(S) = -M * (Y -M(S)), where M * is the adjoint operator of M.

  minSG 1 (L 1 (S)) + G 2 (L 2 (S)) + H(S), One can use the algorithms suggested before and minimization will require the computation of αG * 2 proximity operator, for some given positive real α which is simply given by prox αG * 1 (Z) = Z -Z ,with Z[i, j] = SoftThresh λw j [i] (Z[i, j]).

  min α J (α) s. t. α[l, :] 0 ≤ η l , l = 1 . . . r,(D.3)

J

  (α) s. t. α ∈ E α * . (D.6)The proposed scheme is motivated by the idea of identifying approximately E α * along with the iterative process. One can think of numerous algorithms to solve Problem D.6, all involving the orthogonal projection onto E α * . We build upon the fast proximal splitting algorithm introduced in[START_REF] Beck | [END_REF]. For a vector u ∈ R q we note σ a permutation of 1,q verifying |u[σ(1)]| ≥ • • • ≥ |u[σ(q)]|. For an integer k ≤ q, we define Supp k (u) = {i ∈ 1, q /|u[i]| ≥ |u[σ(k)]|}, (D.7)Finally for a vector α ∈ R r×q , we define the subspaceE k,α = {M ∈ R r×q /Supp k (M[i, :]) = Supp k (α[i, :]), i = 1 . . . r}. (D.8)The proposed scheme is given in Algorithm D.1. f is a positive valued concave increasing function and proj E f (k) ,U k (.) denotes the orthogonal projection onto E f (k) ,U k .Algorithme D.1 Beck-Teboulle proximal gradient algorithm with variable proximity operator 1: Initialization : α 0 = 0 R r×q , β 0 = α 0 , t 0 = 1 res -1 = 0, res 0 = 0, tol, k = 0 2: Minimization 3: tant que k < k max and |(res k -res k-1 )/res k | faire 4:

  res k+1 = J (β k ) 10: k = k + 1 11: Return : α kstop . The solution support size is constraint at step 5 and the size is gradually increased as shown in Fig.D.1. The convergence analysis this scheme is out of the scope of this paper. However, Fig. D.2 suggests that once an index is included in an iterate support, this index is included in all the forthcoming iterates supports. This implies that at each support size's step in Fig. D.1, the algorithm approximately solves a problem of the following form : min α J (α) s. t. α ∈ E, (D.9)

Figure D. 1 -

 1 Figure D.1 -Support size function ; X axis : iteration index k in Algorithm D.1 ; Y axis : f (k) for f (x) = √ x + 1

  ∆ e if i = k + 1,and(definition of the convolution with a zero boundary condition, for a vector v = (v i ) 1≤i≤p , the vector h = v ψ e,a is given by h[j] = p i=1 v i ψ j+k+1-i , (E.3)for j ∈ 1, p and with the convention thatψ j+k+1-i = 0 if j + k + 1 -i < 1 or j + k + 1 -i > p.Combining Eq.E.1, E.2 and E.3, we can writeh[j] = (2 k n=1 a n e ∆ e )v j -i∈[max(1,j-k), min(p,j+k)],i =j 1 |j -i| e ∆ e v i . (E.4) We recall that Ψ e,a (v) = h 2 2 . On the other hand, Ψ e,a (v) = t v e ∆ e v i if j = k + 1 and t v [k + 1-i| e ∆ e v i . (E.5) Thus, t v [k + 1] -h[k + 1] = 0 and for j = k + 1 t v [j] -h[j] = ( e ∆ e v i .(E.6)

  t v [j] -h[j] = ( 1) e ∆ e v 1 -p-j n=k 1 n e ∆ e v n . (E.7)Now, using the inequalities for n > 1, e > 1, we get the following upper bounding :|t v [j] -h[j]| ≤ 1 e -1 [max(|(p -j + 1) 1-e + (j -1) 1-e -(k + 1) 1-e -k 1-e |,|(p -j) 1-e + j 1-e -(k + 1) 1-e -k 1-e |)a + e -1 (j -1) e + k 1-e -(p -j)1-e ]

  Figure G.1 -Point cloud to image mapping.

,

  U 5 = (k 2 -l 2 )1≤k≤N l 1≤l≤Nc , U 6 = (kl)1≤k≤N l 1≤l≤Nc

  Table des matières vii B Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C Numerical Optimal Transport . . . . . . . . . . . . . . . . . . C.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . C.2 General notions . . . . . . . . . . . . . . . . . . . . . . Field-of-view mapping . . . . . . . . . . . . . . . . . . D.3 Barycentric coordinates . . . . . . . . . . . . . . . . . D.4 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . E Numerical results . . . . . . . . . . . . . . . . . . . . . . . . .

	viii		Table des matières
	G Points clouds to image transform
	H Ellipticity parameters directional derivatives
	C.3 Bibliographie	Data representation in the Transport Framework . . .
	D	Transport based PSF field interpolation . . . . . . . . . . . . .
		D.1	Local non-linear dimension reduction . . . . . . . . . .
		D.2
	A Academic activities
	B	Convex analysis
	A Proximal calculus . D Joint super-resolution minimization schemes
	A	Components estimation problem . . . . . . . . . . . . . . . . .
	B	Coefficients estimation . . . . . . . . . . . . . . . . . . . . . .
	E	Notch filter approximation
	F	Multidimensional scaling

E.1 Quality assessment . . . . . . . . . . . . . . . . . . . . E.2 Experiments . . . . . . . . . . . . . . . . . . . . . . . . E.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . F Practical considerations . . . . . . . . . . . . . . . . . . . . . F.1 Parameters . . . . . . . . . . . . . . . . . . . . . . . . F.2 Transportation issues . . . . . . . . . . . . . . . . . . . G Reproducible research . . . . . . . . . . . . . . . . . . . . . . H Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VII Conclusion A Methods and results . . . . . . . . . . . . . . . . . . . . . . . B Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B Convex conjugate . . . . . . . . . . . . . . . . . . . . . . . . . C Multiple frame super-resolution A First guess noise . . . . . . . . . . . . . . . . . . . . . . . . . . B Positivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C Reweighting . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Sommaire A Manifold data model . . . . . . . . . . . . . . . . . . . . . . . 43 A.1 Reminder on differential geometry . . . . . . . . . . .

43 A.2 Manifold learning techniques . . . . . . . . . . . . . . 46 B Optimal transport . . . . . . . . . . . . . . . . . . . . . . . . 52 B.1 Problem formulations . . . . . . . . . . . . . . . . . . 53 B.2 Geometry of optimal transport . . . . . . . . . . . . 55

  that the coordinates slices of the solution of Problem III.13, i.e. the vectors (y i [l]) 1≤i,j≤N , for l ∈ 1, d are solutions of the generalized eigenvalues problem

	Ly = λDy,	(III.15)

  Sommaire A PSF field interpolation state-of-the-art . . . . . . . . . . . . . 108 B Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 C Numerical Optimal Transport . . . . . . . . . . . . . . . . . 109 C.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . 110 C.2 General notions . . . . . . . . . . . . . . . . . . . . . 110 C.3 Data representation in the Transport Framework . . 113

D Transport based PSF field interpolation . . . . . . . . . . . . 114 D.1 Local non-linear dimension reduction . . . . . . . . . 115 D.2 Field-of-view mapping . . . . . . . . . . . . . . . . . 116 D.3 Barycentric coordinates . . . . . . . . . . . . . . . . . 117 D.4 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 118 E Numerical results . . . . . . . . . . . . . . . . . . . . . . . . 119 E.1 Quality assessment . . . . . . . . . . . . . . . . . . . 119 E.2 Experiments . . . . . . . . . . . . . . . . . . . . . . . 120 E.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 122 F Practical considerations . . . . . . . . . . . . . . . . . . . . . 125 F.1 Parameters . . . . . . . . . . . . . . . . . . . . . . . . 125 F.2 Transportation issues . . . . . . . . . . . . . . . . . . 128 G Reproducible research . . . . . . . . . . . . . . . . . . . . . . 137 H Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

  a 1 + a 2 t c t -e 1 (X i + tP j )

					d 1 + d 2 t c t	(H.3)
	de 2 (X i + tP j ) dt	=	b 1 + b 2 t c t	-e 2 (X i + tP j )	d 1 + d 2 t c t	,	(H.4)

  Georges Oppenheim, editeurs, Wavelets and Statistics, volume 103 of Lecture Notes in Statistics, pages 55-81. Springer New York, 1995. (cited on page 75) [Burkard et al. 2009] Rainer Burkard, Mauro Dell'Amico et Silvano Martello. Assignment problems. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2009. (cited on page 112) [Cabral et al. 2013] Ricardo Cabral, Fernando De La Torre, Joao P Costeira et Alexandre Bernardino. Unifying nuclear norm and bilinear factorization approaches for low-rank matrix decomposition. In Proceedings of the IEEE International Conference on Computer Vision, pages 2488-2495, 2013. (cited on page 40) [Candès & Donoho 2004] Emmanuel J Candès et David L Donoho. New tight frames of curvelets and optimal representations of objects with piecewise C2 singularities. Communications on pure and applied mathematics, vol. 57, no. 2, pages 219-266, 2004. (cited on page 36) [Candès & Plan 2010] Emmanuel J Candès et Yaniv Plan. Matrix completion with noise. Proceedings of the IEEE, vol. 98, no. 6, pages 925-936, 2010. (cited on page 40) [Candès & Recht 2009] Emmanuel J Candès et Benjamin Recht. Exact matrix completion via convex optimization. Foundations of Computational mathematics, vol. 9, no. 6, pages 717-772, 2009. (cited on page 40) [Candès & Romberg 2007] Emmanuel Candès et Justin Romberg. Sparsity and incoherence in compressive sampling. Inverse problems, vol. 23, no. 3, page 969, 2007. (cited on page 37) [Candès et al. 2006] Emmanuel J Candès, Justin K Romberg et Terence Tao. Stable signal recovery from incomplete and inaccurate measurements. Communications on pure and applied mathematics, vol. 59, no. 8, pages 1207-1223, 2006. (cited on page 34) [Candès et al. 2008] Emmanuel J Candès, Michael B Wakin et Stephen P Boyd. Enhancing sparsity by reweighted 1 minimization. Journal of Fourier analysis and applications, vol. 14, no. 5-6, pages 877-905, 2008. (cited on pages 39, 66, 67 et 93) [Candès et al. 2011] Emmanuel J Candès, Xiaodong Li, Yi Ma et John Wright. Robust principal component analysis ? Journal of the ACM (JACM), vol. 58, no. 3, page 11, 2011. (cited on page 40) [Cartis & Thompson 2015] C. Cartis et A. Thompson. A New and Improved Quantitative Recovery Analysis for Iterative Hard Thresholding

http://terapix.iap.fr/cplt/T0007/doc/T0007-doc.html

http://www.astro.washington.edu/docs/idl/htmlhelp/slibrary21.html

Academic activities

Journal papers -F. M. Ngolè Mboula, J.-L. Starck, S. Ronayette,K. Okomura, J. Amiaux.

Super-resolution method using sparse regularization for point spread function recovery, Astronomy and Astrophysics, 575 : A86, 2015.

-F. M. Ngolè Mboula, J.-L. Starck, K. Okomura, J. Amiaux, P. Hudelot.

Constraint matrix factorization for space variant PSFs field restoration, Inverse Problems, 32 : 12, 2016.

-F. M. Ngolè Mboula, J.-L. Starck. PSFs field learning based on Optimal Transport distances, Submitted.

-S. Farrens, F. M. Ngolè Mboula, J.-L. Starck. Space variant deconvolution of galaxies survey images, Astronomy and Astrophysics.

-Great3 Challenge participants. GREAT3 results-I. Systematic errors in shear estimation and the impact of real galaxy morphology, Monthly Notices of the Royal Astronomical Society, 450(3) : 2963-3007, 2015.

Conferences and workshops

-International Traveling Workshop on Interactions between Sparse models and Technology, Namur, Belgium, August 2014 (poster).

-International Workshop on Cosmology and Sparsity, Nice, September 2014 (talk).

-RICAM Special Semester on New Trends in Calculus of Variations : Optimal Transport in the Applied Sciences, Linz, December 2014 (talk).

-Data Sience in Particle Physics, Astrophysics and Cosmology workshop, Paris, April 2015 (talk).

-Missing data in physics, Nice, May 2015 (talk).

-Euclid weak lensing meeting, London, November 2015 (talk).

-International Traveling Workshop on Interactions between Sparse models and Technology, Aalborg, Denmark, August 2016 (talk).

Annexe B. Convex analysis

Chapter IV Section B.1, and for x ∈ R d 2 p 2 , we have,

-Another popular example of proximity operator is the one associated with F = λ . 1 with λ ∈ R :

where p is the dimension of H, α i the components of alpha in the basis associated with . 1 and (.) + = max(., 0).

-For λ = (λ i ) 1≤i≤p , the proximity operator associated with the weighted

The hard thresholding operator defined as

is often used instead in practice.

Properties :

-Now we suppose that H = R p , and we want to solve

where

Many problems in signal and image processing may be formulated this way, where F 1 would be the data attachment function Annexe F

Multidimensional scaling

We consider a family of vectors (r i ) 1≤i≤p in R q ; we set R = [r 1 , . . . , r p ]. We assume that p i=1 r i = 0 R q . We define the distances matrix D = ( r i -r j 2 2 ) 1≤i,j≤p and the norms vector v = ( r 1 2 2 , . . . , r p 2 2 ) T . We want to retrieve (r i ) 1≤i≤p from D. D can be rewritten as :

We consider the centering matrix C introduced in Chapter VI, Section D.1. We have

Points clouds to image transform

The step 11 in the Algorithm VI.1 in Chapter VI.1 consists in converting estimated Wasserstein barycenters which are points clouds into images which are the actual PSFs estimates. Let consider an estimated Wasserstein X bar ∈ R 3×N . By construction, the first dimension is related to pixels intensities and the second and third are related to pixels positions (see Eq. VI.12). We assume that the PSFs images have N l lines and N c columns so that N = N l N c . Then the most simple way of calculating an image from X bar is the following :

1. Initialize an N l × N c image PSF bar with all pixels values set to zeros.

FOR ALL

However, this generally results in sharp pixels intensities variations and therefore visual artifacts. We circumvent this by involving the four nearest pixels in step 2 in the procedure above, rather than the nearest one solely :

1. Initialize an N l × N c image PSF bar with all pixels values set to zeros. Résumé : L'un des challenges majeurs de la cosmologie moderne réside en la nature même de la matière et de l'énergie noire. La matière noire peut être directement tracée à travers son effet gravitationnel sur les formes des galaxies. Abstract : One of the biggest challenges of modern cosmology is to gain a more precise knowledge of the dark energy and the dark matter nature. Fortunately, the dark matter can be traced directly through its gravitational effect on galaxies shapes. The European Spatial Agency Euclid mission will precisely provide data for such a purpose. A critical step is analyzing these data will be to accurately model the instrument Point Spread Function (PSF), which the focus of this thesis.

FOR ALL

This is illustrated in

We developed non parametric methods to reliably estimate the PSFs across an instrument field-ofview, based on unresolved stars images and accounting for noise, undersampling and PSFs spatial variability. At the core of these contributions, modern mathematical tools and concepts such as sparsity.

An important extension of this work will be to account for the PSFs wavelength dependency.

Université Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France