Gabriele Facciolo
email: gabriele.facciolo@upf.edu

Andrés És Almansa
email: andres.almansa@telecom-paristech.fr

Jean-François Aujol
email: jean-francois.aujol@cmla.ens-cachan.fr

Vicent Caselles
email: vicent.caselles@upf.edu

Jean-Franc ¸ois

Irregular to regular sampling, denoising and deconvolution

Keywords: Image restoration, Total Variation, variational methods, satellite images AMS subject classifications. 68U10, 65K10, 65J20, 94A08

come

Introduction.

A general image acquisition system may be modeled by the following image formation model

z(ξ k ) = (h * u)(ξ k ) + n ξ k , ξ k ∈ Ξ, (1.1)
where Ξ = {ξ k } N 2 k=1 ⊆ R 2 is a finite set of regular or irregular samples, u : R 2 → R is the ideal undistorted image, h : R 2 → R is a blurring kernel whose Fourier spectrum ĥ has most of its energy concentrated in the spectral support of u, z is the observed sampled image which is represented as a function z : Ξ → R, and n ξ k is, as usual, a white Gaussian noise with zero mean and standard deviation σ.

Reconstructing a signal u : R 2 → R over an infinite support from a finite set of samples z(ξ k ) is not possible without imposing restrictions. As in most works, in order to simplify this problem, we shall assume that the functions h and u are periodic of period N in each direction. Let us denote by Ω N the interval [0, N [ 2 and assume that h, u are functions defined in Ω N . To fix ideas, we assume that h, u ∈ L 2 (Ω N ), so that h * u is a continuous function in Ω N [START_REF] Gasquet | Analyse de Fourier et Applications[END_REF] (which may be extended to a continuous periodic function in R 2 ). Then the samples (h * u)(ξ k ), ξ k ∈ Ξ, are well defined.

We shall concentrate on the particular case of perturbed sampling and we shall assume that Ξ is a set of N 2 samples, which take the particular form

Ξ = Z 2 ∩ Ω N + ε(Z 2 ∩ Ω N ), (1.2)
where ε : R 2 → R 2 is a "smooth and small" perturbation function in the sense that supp ε ⊆ [-N Tε , N Tε ] 2 for some period T ε > 2 corresponding to the maximum vibration frequency and the mean amplitude of the perturbation ([ ΩN |ε(x)| 2 dx/ ΩN dx]

1 2 ) is small with respect to 1 pixel (we refer to Section 2 for a model (2.1) of this perturbation and also for a general overview of irregular sampling aspects).

Even knowing the exact sampling geometry Ξ, the blurring kernel h and the statistics of the noise n, the problem of recovering u from z is very ill-posed in a variety of situations that arise in real-world image acquisition systems. The main difficulty comes from the fact that perfect reconstruction from irregular samples can only be ensured either with a sampling rate, which is much larger than in the regular case (relative to the bandwidth of h), or similar to the regular critical sampling rate, but under very constrained sampling geometries that rarely appear in practice.

(See [START_REF] Almansa | Echantillonnage, Interpolation et Détection[END_REF]Chapter 2] for a discussion on this problem). In any case, even if perfect reconstruction is sometimes possible in theory, with a relatively small and realistic sampling rate, it is still very ill posed, which limits its practical interest whenever noisy measurements are involved. A common and unavoidable strategy to solve this ill-conditioning is regularization. And the typical constrained formulation of a regularization method [START_REF] Tikhonov | Solutions of Ill-Posed Problems[END_REF] consists in choosing between all possible solutions of (1.1) the one that minimizes the functional J(u). Modeling the acquisition system (1.1) through constrained minimization eliminates the regularization parameter in classical variational or bayesian formulations, and due to the statistical nature of the noise information, the constraint is expressed as an upper estimate of the noise variance σ 2 . The constrained restoration problem becomes min u J(u), subject to

ξ k ∈Ξ |(h * u)(ξ k ) -z(ξ k )| 2 ≤ N 2 σ 2 ,
where the regularizer J(u) embodies our a-priori knowledge of the image, specifying its smoothness properties. The use of the Dirichlet integral ΩN |Du| 2 dx is not satisfactory, mainly due to the unability of the previous functional to resolve discontinuities (edges) and oscillatory textured patterns, the information corresponding to high frequencies of z being attenuated by it. Indeed, functions in W 1,2 (Ω) (i.e., functions u ∈ L 2 (Ω) such that Du ∈ L 2 (Ω)) cannot have discontinuities along rectifiable curves. These observations motivated the introduction of Total Variation (T V (u) = ΩN |Du|) in image restoration problems by L. Rudin, S. Osher and E. Fatemi in their work [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. The a priori hypothesis is that functions of bounded variation (the BV model) ( [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF][START_REF] Evans | Measure Theory and Fine Properties of Functions[END_REF]) are a reasonable functional model for many problems in image processing, in particular, for restoration problems ( [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]). Typically, functions of bounded variation have discontinuities along rectifiable curves, being continuous in the measure theoretic sense away from discontinuities. The discontinuities could be identified with edges. The ability of total variation regularization to recover edges is one of the main features which advocates for the use of this model (its ability to describe textures is less clear, some textures can be recovered, but up to a certain scale of oscillation). We refer to [START_REF] Evans | Measure Theory and Fine Properties of Functions[END_REF] for the definition of functions of bounded variation and their basic properties.

We shall explore in this paper a family of regularizers that takes into account the spectral decay of the Fourier coefficients in the class of images we are looking for. In the case of satellite images, this spectral behavior can be estimated by statistical measures of the decay of Fourier coefficients. The general class of regularizers we consider is

J A (u) = ΩN |A(D)u| (1.3)
where A(D)u is defined by the coefficients of its Fourier series F(A(D)u)(ω) = A(iω)û(ω), ω ∈ Z 2 . Note that J A (u) < ∞ imposes a frequency penalization according to the profile A(iω). In practice we choose A(iω) so that |A(iω)| ∼ | 2π N ω| p for large |ω|, 1 ≤ p ≤ 2. This is in consonance with the approach of Gröchenig and Strohmer [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF] that proposes to incorporate an a-priori decay in the restoration process (see Section 2).

In [START_REF] Almansa | A TV based restoration model with local constraints[END_REF] (see also [START_REF] Almansa | Restoration and zoom of irregularly sampled, blurred and noisy images by accurate total variation minimization with local constraints[END_REF][START_REF] Bertalmío | TV based image restoration with local constraints[END_REF][START_REF] Rudin | Total Variation based Image Restoration with Free Local Constraints[END_REF]), the authors proposed a restoration model using a local estimate of the noise variance. The local formulation contributes to reducing the unnaturally-looking aspect of images obtained from global Total Variation based minimization, thus improving the recovery of textures. Following the mentioned proposal we replace the constraint

ξ k ∈Ξ |(h * u)(ξ k ) -z(ξ k )| 2 ≤ N 2 σ 2 , by G * |∆ Ξ (h * u) -z| 2 (ξ k ) ≤ σ 2 , ∀ξ k ∈ Ξ, (1.4)
where the sampling operator ∆ Ξ : C(R 2 ) → ℓ 2 (Ξ) is given by ∆ Ξ (v) = {v(ξ k )} N 2 k=1 and G is a discrete convolution kernel such that G(ξ) > 0 for all ξ ∈ Ξ and k G(ξ k ) = 1.

Combining the two ideas described above, the use of a regularizer that takes into account the spectral decay of images in a certain class (1.3), and the incorporation of the image acquisition model as a set of local constraints (1.4), we propose the following constrained variational model for restoring u

min u Ω |A(D)u|, subject to G * |∆ Ξ (h * u) -z| 2 (ξ k ) ≤ σ 2 ∀ξ k ∈ Ξ.
(1.5)

The constrained formulation (1.5) can be solved using the unconstrained formulation

min u max (λ k )≥0 Ω |A(D)u| + 1 2 ξ k ∈Ξ λ k G * |∆ Ξ (h * u) -z| 2 (ξ k ) -σ 2 (1.6)
where λ k ≥ 0 is a Lagrange multiplier that has to be chosen so that the constraints (1.4) are satisfied. Let us say explicitly that both the blurring kernel h and the sampling grid Ξ (alternatively the grid perturbation function ε) are assumed to be known exactly, and that the only thing known about the noise n ξ k is that it is a white Gaussian noise with zero mean and known variance σ 2 . Several methods exist to estimate all these parameters [START_REF] Jalobeanu | Estimation of blur and noise parameters in remote sensing[END_REF] for a given acquisition device and we shall not address this question here. The case of recovering an irregularly sampled image on a regular sampling grid was already considered by the second author in [START_REF] Almansa | Echantillonnage, Interpolation et Détection[END_REF], but the blurring kernel h was assumed to be an ideal window (with Nyquist frequency cutoff), i.e., ĥ = χ [-1/2,1/2] 2 . Different numerical algorithms were tested in the case where the sampling set is perturbed according to (1.2) and they worked relatively well only within a low-frequency spectral region R ⊆ [-α, α] 2 , where α ≈ 1 2 -1/T ε . When attempting to recover û in the high frequency band [-1/2, 1/2] 2 \R serious theoretical and numerical problems appeared and, actually, restoration errors were most important there. Subsequently, the restoration problem (1.5) was studied in [START_REF] Almansa | Restoration and zoom of irregularly sampled, blurred and noisy images by accurate total variation minimization with local constraints[END_REF] when J(u) is the total variation and the image acquisition model was incorporated as a set of local constraints on a partition of the image obtained as a result of a segmentation. The use of local constraints (1.4) was advocated in [START_REF] Almansa | A TV based restoration model with local constraints[END_REF] and we also adopt this technique here.

Let us finally mention that many numerical algorithms have been proposed to minimize total variation (or similar models) subject to a global constraint as in (1.5) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Geman | Constrained Image Restoration and Recovery of Discontinuities[END_REF][START_REF] Blanc-Féraud | Nonlinear image processing: modeling and fast algorithm for regularization with edge detection[END_REF][START_REF] Vogel | Iterative methods for total variation denoising[END_REF][START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF][START_REF] Chan | Continuation method for total variation denoising problems[END_REF][START_REF] Chan | On the convergence of the Lagged Diffusivity Fixed Point Method in total variation image restoration[END_REF][START_REF] Charbonnier | Deterministic edgepreserving regularization in computed imaging[END_REF][START_REF] Durand | Image Deblurring, Spectrum Interpolation and Application to Satellite Imaging[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Scherzer | Variational Methods in Imaging[END_REF]. Imposing local constraints in a partition of the image was proposed in [START_REF] Rudin | Total Variation based Image Restoration with Free Local Constraints[END_REF] and further developed in [START_REF] Bertalmío | TV based image restoration with local constraints[END_REF][START_REF] Almansa | Restoration and zoom of irregularly sampled, blurred and noisy images by accurate total variation minimization with local constraints[END_REF][START_REF] Almansa | A TV based restoration model with local constraints[END_REF]. In [START_REF] Lintner | Solving a variational image restoration model which involves L ∞ constraints[END_REF] the authors combined total variation minimization with a set of constraints of type | h * uz, ψ | ≤ τ where ψ varies along an orthonormal basis of wavelets (or a family of them) and τ > 0. The aim was also to construct an algorithm which preserves textures and has good denoising properties. As we will do here, these constraints were incorporated using Uzawa's algorithm. In [START_REF] Gilboa | Texture Preserving Variational Denoising Using an Adaptive Fidelity Term[END_REF], the authors proposed to minimize total variation subject to a family of local constraints which control the local variance of the oscillatory part of the signal. The constraints are introduced via Lagrange multipliers with an approach similar to the one used in [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. This amounts to adding a spatially varying fidelity term that locally controls the extent of denoising over image regions depending on their content. Besides the fact that we use Uzawa's algorithm and we try to address the problem of deconvolution and denoising of irregularly sampled images, the work [START_REF] Gilboa | Texture Preserving Variational Denoising Using an Adaptive Fidelity Term[END_REF] is quite similar to our approach. Finally, in [START_REF] Scherzer | Variational Methods in Imaging[END_REF], the authors proposed a non-convex data attachment term with a larger weight (depending on the inverse of the modulus of the gradient) in flat areas than in textured ones. This goes in the same sense as our model, that is, to be able to denoise flat regions of the image while keeping the oscillations in textured areas.

Let us summarize the contributions of the paper. Our goal is to propose an algorithm for image denoising and deconvolution. The specificity of the considered problem comes from the fact that the images we deal with are irregularly sampled. On the one hand, we carry out a thorough and accurate modelization of the problem based on local constraints (1.4), and propose a family of frequency adaptive regularization functionals (1.3) that incorporates a priori image smoothness model. On the other hand, we prove the convergence of our algorithm. Lastly, the study of the local constraint model combined with standard error propagation allows to derive useful stopping criteria for the algorithm.

Let us finally explain the plan of the paper. In Section 2 we introduce the problem of irregular to regular sampling and we discuss the ACT algorithm of Gröchenig and Strohmer [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF]. In Section 3 we introduce our frequency adaptive variational restoration model with local constraints in the continuous setting. The corresponding discrete model is discussed in Section 4 along with the adaptation of a computational improvement introduced by L. Moisan in [START_REF] Moisan | How to discretize the Total Variation of an image?[END_REF]. In Section 5 we study the existence, uniqueness and numerical approximation to the model introduced in the previous Section. This study is completed in Section 6 where we describe a Quasi-Newton algorithm for the solution of the Euler Lagrange equation corresponding to the energy in (1.6) for fixed values of the Lagrange multipliers (λ k ). In Section 7 we propose a practical stopping condition for the restoration algorithm for the local constraint model. In Section 8 we display experiments concerning restoration and zooming of irregularly sampled images. Section 9 summarizes the main conclusions of this work.

1.1. Preliminaries and notations. Let us introduce some notations.

For any function u ∈ L 2 (Ω N ) (assuming periodicity of period N in each direction) we denote its Fourier coefficients as

û(p, q) = 1 N 2 ΩN u(x, y)e -2πi (px+qy) N dxdy for (p, q) ∈ Z 2 .
As in [START_REF] Malgouyres | Edge direction preserving image zooming: a mathematical and numerical analysis[END_REF], our plan is to compute a band limited approximation to the solution of the restoration problem. To do so, assume for simplicity that M is an even number and define

B M := {u ∈ L 2 (Ω N ) : û is supported in I M } where I M := {- M 2 + 1, . . . , M 2 } 2 .
We notice that B M is a finite dimensional vector space of dimension M 2 which can be identified with R M 2 by mapping u ∈ B M to the vector

u = u(r N M , l N M ) M-1 r,l=0 . Moreover, if u ∈ B M we may write u(x, y) = -M 2 <p,q≤ M 2 û(p, q)e 2πi (px+qy) N .
where

û(p, q) = 1 M 2 0≤r,l<M u r N M , l N M e -2πi (pr+ql) M , - M 2 < p, q ≤ M 2 .
Then the values u r N M , l N M , 0 ≤ r, l < M , can be recovered as the discrete inverse Fourier transform of û(p, q). Hence u ∈ B M can also be identified with the vector of Fourier coefficients û ∈ I C M 2 . Note that we shall mainly study here the critical sampling case M = N , however we will keep two different symbols M for the bandwidth and N for the domain size in order to simplify the exposition of certain parts of the algorithm.

We intend to solve the restoration problem in the class of band-limited functions B M . Later on we will comment on this choice. We will also use the operator notation for the Fourier transform that applied to the function u returns the vector of its Fourier coefficients: û = Fu , conversely F * û = u denotes the inverse transform, then

F * F = Id.
2. Irregular to regular sampling. Opposed to digital photographs, satellite images are generally not acquired by a squared array of sensors but by a sweeping bar of sensors known as TDI (Time Delay Integrator) [START_REF] Schott | Remote Sensing: The Image Chain Approach[END_REF]. This acquisition geometry called push-broom is widely applied in aerospace imaging applications and, nowadays, it provides the highest resolution in earth imaging applications. As a consequence of this progressive acquisition mode, the micro-vibrations of the satellite together with irregularities in sensors position result in perturbed sampling sets. In most cases, the knowledge of certain vibration modes and the analysis of acquired images help to estimate, very accurately, the perturbations in the sampling grid, which can be modeled [START_REF] Almansa | Echantillonnage, Interpolation et Détection[END_REF] by

ε(x) = q k=1 a k (x)cos(2π ω k N , x + φ k ), x ∈ R 2 , (2.1)
for some q ≥ 1, where a k (x) are smooth modulation functions and the vibration frequencies ω k are an order of magnitude (or even more) below the Nyquist frequency of the sampling rate. The bound on the modulation functions is inversely proportional to ω k and the number of vibration modes is small. This results in smooth and small perturbations with |ε(x)| no larger than a few pixels, and perturbation slope |∇ε(x)| no larger than about one tenth of a pixel per pixel. As a consequence these perturbations are hardly noticeable and we should talk of perturbed sampling rather than irregular sampling in those cases. Even if the image distortion is not evident from a geometrical point of view it is very important to correct the perturbations in image registration applications where a sub-pixel accuracy is necessary. In order to be less dependent on a particular physical instrument, in our experiments we used a simplified version of this model which still captures its main characteristics, namely the perturbation function ε = (ε 1 , ε 2 ) is simulated as a discrete colored noise, i.e. for ω ∈ Z 2 we define

ε i (ω) ∼ N (0, σ2 ) if |ω| ≤ N/T ǫ = 0 otherwise (2.2)
where σ is chosen in such a way that the standard deviation of ε i (x) is A for i ∈ {1, 2}. This gives σ = ATǫ 2 (we have taken the Fourier transform as an isometry). Thus the behavior of the perturbation is characterized by the two parameters "amplitude" A and maximal vibration frequency N/T ε (or "minimal vibration period" T ε ). The precise values of A and T ε used in our experiments will be specified in the experiments Section.

There are many works in the literature dealing with the irregular to regular sampling problem. Let us mention that, according to Kadec's theorem [START_REF] Kadec | The exact value of the Paley-Wiener constant[END_REF], we have a perfect recovery of the signal if we consider a perturbed sampling with small perturbations |ε(x)| ≤ 0.11. Recall also that Beurling-Landau's theorem [START_REF] Landau | Necessary density conditions for sampling and interpolation of certain entire fuctions[END_REF], ensures perfect reconstruction of a function from its samples for arbitrary stable sampling sets [START_REF] Landau | Necessary density conditions for sampling and interpolation of certain entire fuctions[END_REF], but it requires the (lower) sampling density to be larger than 1. These conditions are very restrictive and do not hold true for most of the image restoration applications. For a comparison between several iterative methods we refer the reader to [START_REF] Azmi | Comparison between several iterative methods of recovering signals from nonuniformly spaced samples[END_REF][START_REF] Feichtinger | Iterative algorithms in irregular sampling: A first comparison of methods[END_REF][START_REF] Almansa | Echantillonnage, Interpolation et Détection[END_REF].

2.1. The ACT algorithm. One of the best performing reconstruction methods available for irregular to regular sampling is the ACT algorithm (for Adaptive-weights Conjugate-gradient on Toeplitz-matrix) introduced by Gröchenig et. al. in [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF]. This method represents a discrete image u as a trigonometric polynomial of order M/2 in each variable (for simplicity of notation we shall assume that M is an even number) so that the interpolation at the sampling points Ξ

= {ξ k } N 2 k=1 ⊆ R 2 becomes u(ξ k ) = t∈{-M 2 +1,..., M 2 } 2 a t e 2πi M t,ξ k , k ∈ {1, . . . , N 2 }. (2.3)
Thus, if z represents the irregularly sampled data we may write [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF] z = Sa, where S = ((s kt )),

s kt = e 2πi M t,ξ k , (2.4)
i.e. S is the Vandermonde matrix associated to the trigonometric polynomial in (2.3). Note that S maps a ∈ ℓ 2 ({-M 2 + 1, . . . , M 2 } 2 ) to z = {z(ξ k )} ∈ ℓ 2 (Ξ) as given in (2.3). The bandwidth M of the trigonometric polynomial is chosen to be M ≤ N , so the system (2.4) is expected to be determined or over-determined.

Following [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF], the ACT algorithm recovers the coefficients a of the trigonometric polynomial by solving the least squares problem arg min

a √ W (Sa -z) 2 , (2.5)
where the matrix W = diag({w k } k=1..N 2 ) assigns weights that are inversely proportional to the sampling density at ξ k :

w k = area(V k ) where V k := ({x : |x -ξ k | < |x -ξ j |, ∀j = k}) . (2.6)
If we interpret the discrete data z(ξ k ) as a piecewise constant function k z(ξ k )χ V k , then the weights w k guarantee the isometry between the irregular sampling on the grid Ξ and its function representation, thus compensating the local variations in the sampling density. Moreover, by using the weights W , Gröchenig and Strohmer provide an explicit estimate for the rate of convergence of the ACT algorithm [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF].

The system of normal equations associated to (2.5) is

S * W Sa = S * W z , (2.7)
where S * denotes the adjoint matrix of S (a notation that will be used throughout this paper). Observe that the M 2 × M 2 matrix S * W S has a Toeplitz structure [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF], and thus, S * W Sa is efficiently computed in O(M 2 log 2 (M 2 )) time using Fourier methods. Moreover the entries of T := S * W S and b := S * W z can be approximated using the NFFT [START_REF] Potts | Fast Fourier transforms for nonequispaced data: A tutorial[END_REF] in O(M 2 log 2 (M 2 )) time each [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF]. Finally, (2.7) is solved using a conjugate gradient (CG) method. The following algorithm summarizes the method.

Algorithm I: ACT algorithm for a fixed bandwidth M Requires: N 2 irregular samples in vector z.

Ensures: M 2 regular samples in vector u.

1. Compute T = S * W S and b = S * W z using the NFFT.

2. Solve T a = b using conjugate gradients.

3. Compute the regular samples u(i N M , j N M ) for (i, j) ∈ {0, . . . , M -1} 2 by applying the inverse FFT to a.

Let us note that in the more realistic cases where T is not invertible or illconditioned, the CG solver acts as a regularizer and chooses the minimum norm solution a among those that satisfy (2.7). This is a constrained variational formulation that can be written using a Lagrange multiplier λ > 0 as the unconstrained minimization problem

min a a 2 + λ T a -b 2 .
This formulation also applies to the following two variants of the ACT algorithm [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF] that incorporate an a-priori spectral decay rate |û(ω)| ≤ Lφ(ω), for some L > 0, of the image class to be restored (when available). For satellite images this estimation has been performed by Almansa in [START_REF] Almansa | Echantillonnage, Interpolation et Détection[END_REF] and it is given by φ(ω) = (1 + | 2π N ω|) -p for some value of p near 1.6.

• The first proposed variant solves CT a = Cb (instead of (2.7)) using the CG algorithm, where

C = diag({φ(ω)} ω∈{-M 2 +1,..., M 2 } 2
). Notice that we can write the problem as

(ACT D ) min a a 2 + λ C(T a -b) 2 , (2.8)
for some λ > 0. Since the weighting matrix C is applied to the residuals T ab, then it affects the relative cost of the errors by shifting the cost from high frequencies towards the lower ones. This shift is reflected in CG search directions, altering the intermediate solutions so that they will fit the low frequency before the higher ones. • The second ACT variant solves T d = b where a = Cd, using the CG algorithm. Re-writing it as an optimization problem, observe that the weights appear now in the regularity term

(ACT R ) min a C -1 a 2 + λ T a -b 2 .
(2.9)

The spectral weights C -1 are now penalizing the apparition of higher frequencies in the solution a and not in the residual. In either case, if T is invertible and the CG algorithm converged, then the solution of both variants coincides with the solution of (2.7). But the CG iteration is truncated before its convergence mainly due to the ill-conditioning of the operator T . So, the solutions obtained by the above methods differ because of the different search directions. As it can be observed experimentally, incorporating the spectral decay indeed reduces the restoration errors, specially when applied to the regularity term in (2.9) (see Table 8.1 in Section 8). In that case, it amounts to finding a solution in a class of functions with a particular spectral decay.

Remark 1. (Global constraint & stopping condition) In [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF] the authors also proposed to extend the ACT algorithm in order to consider the presence of Gaussian noise n (with standard deviation σ) in the image formation model: Sa + n = z. This extension is implemented as a stopping condition for the CG algorithm (Step 2 of Algorithm I)

stop CG if Sa -z 2 ≤ τ N 2 σ 2 with τ ≃ 1
that is designed to avoid the over-fitting of the solution and hence the amplification of the noise. Thus, the ACT with residual-based stop conditions can be seen as a numerical approximation to the minimization of Ca 2 subject to the constraints above. Remark 2. (Fourier vs. spline models) Notice that the ACT Algorithm is based on the underlying assumption that the data can be represented by a trigonometric polynomial. Other interpolation models like the B-Spline have been used in the literature [START_REF] Arigovindan | Variational image reconstruction from arbitrarily spaced samples: a fast multiresolution spline solution[END_REF]. In this work we will restrict ourselves to trigonometric polynomials mainly because convolutions are more easily modeled in this setting, but we intend to explore the use of B-splines in the future.

A frequency adaptive restoration model with local constraints.

In this Section we sketch main elements of the restoration model adopted in this paper. To simplify the presentation the model is first defined in a continuous setting, and the discretization details are deferred to Section 4. The proposed model is inspired from the ideas presented in the previous Section on the ACT, as well as on a recently proposed ACT+TV extension (next Subsection), and the use of local constraints (Subsection 3.3) as in [START_REF] Almansa | A TV based restoration model with local constraints[END_REF].

3.1. The ACT + TV extension. In [START_REF] Almansa | Restoration and zoom of irregularly sampled, blurred and noisy images by accurate total variation minimization with local constraints[END_REF] the authors proposed to combine the ACT algorithm written as (2.5) with total variation regularization, i.e., min

u ΩN |Du|, subject to √ W (S û -z) 2 ≤ N 2 σ 2 . (3.1)
For convenience, let us refer to this model as ACT+TV. As reported in [START_REF] Almansa | Restoration and zoom of irregularly sampled, blurred and noisy images by accurate total variation minimization with local constraints[END_REF] (see also Table 8.1 in Section 8) we observe an improvement of ACT+TV with respect to the original ACT algorithm in terms of MSE error, this improvement is mainly attributed to the edge preserving ability of the TV regularizer.

Frequency Adaptive Regularization (FAR).

Inspired by (2.9) we propose to integrate the spectral weight priors given by the matrix C into the ACT+TV formulation.

From the ACT+TV formulation we keep: (i) the use of the spatial L 1 -norm of a (pseudo)differential operator as a regularizer (instead of the L 2 -norm C -1 a in (2.9)), because of the better edge-preserving capabilities of the L 1 -norm; and (ii) the formulation as an optimization problem under constraints, which enables us to choose automatically the regularization parameter λ in (2.9).

From the ACT R we keep the idea of using a frequency adapted regularizer. For that we define the function

ω ∈ Z 2 → A(iω) ∈ I C × I C. We assume that A(0) = 0, A(iω) = 0 ∀ω = 0, and |A(iω)| ≤ L 1 + 2π N ω κ , ∀ω ∈ Z 2 , (3.2) for some L > 0, κ ≥ 0. If u ∈ C ∞ (Ω N
) can be extended as a smooth and periodic function to R 2 , we define A(D)u by its Fourier coefficients

A(D)u(ω) = A(iω)û(ω) ω ∈ Z 2 .
We define the regularizer functional

J A (u) = ΩN |A(D)u| = ΩN |(A(iω)û(ω)) ∨ |.
Notice that J A (u) can be defined for any u ∈ L 2 (Ω N ) such that A(D)u is a Radon measure. In practice this will not be a problem since we only consider band limited functions.

The total variation J(u) = ΩN |Du| corresponds to the choice A(iω) = 2π N iω. Recall that functions with finite total variation are a good model for image restoration since they permit to recover the discontinuities of the image. But, in practice, digital images may exhibit a stronger decay in its Fourier coefficients than | 2π N ω| -1 and other functional models can be better suited. The operator A(D) permits us to penalize the frequencies according to the profile A(iω). Notice that

|A(iω)û(ω)| = | A(D)u(ω)| ≤ ΩN |A(D)u| =: L for all ω ∈ Z 2 , hence |û(ω)| ≤ L |A(iω)| . If |A(iω)| ∼ | 2π N ω| κ for large |ω|, then |û(ω)| decreases as | 2π N ω| -κ .
In this way we can favor a specific decay rate of the Fourier coefficients of u as a prior.

Remark 3. Notice that the frequency adaptive regularizer of ACT R has the form

C -1 û = φ -1 (ω)û(ω) where φ(ω) = (1 + |2πω/N |) -k , k ≥ 0. In practice, we consider functions A(iω) = 2πiω N ψ(ω) where ψ(ω) = (1 + |2πω/N |) p , p ≥ 0 (p = k -1
gives the same decay). Note that we can write A(D)u = DΨu, where Ψu = ψû, and A(D)u is the derivative of a filtered version of u. Since |A(D)u| 2 = ω |A(iω)û(ω)| 2 , if we use ΩN |A(D)u| 2 as regularizer, then the phase information disappears. This is not the case when we use the L 1 -norm of A(D)u. Besides keeping the structure of a regularization operator, when using the L 1 -norm, this phase is relevant in practice (we obtain bad results if we dismiss it).

Local constraints.

Following [START_REF] Almansa | A TV based restoration model with local constraints[END_REF] we will interpret the acquisition model as a set of local constraints (1.4) and we will search for a solution of the restoration problem that minimizes the frequency adaptive regularizer J A (u) with a data-fitting term derived from the acquisition model and consisting in a set of local constraints (as in (1.5)). For convenience we will describe these constraints directly in the discrete model (see the next Section).

4.

A discrete regularization model. Let us stress here the fact that the regularization functional J A defined in the previous Section is adapted to the restoration of functions with infinite resolution while its numerical approximation restricts the solution to be in a finite dimensional space. We are going to adopt here the following practical point of view. Since our data consists of a finite set of samples, we are going to reconstruct a sampled version of the image and therefore we work in a finite dimensional space. This reflects the fact that digital images have finite resolution. Usually, the restored image is modeled as a piecewise constant function (the values given on the set of pixels), but we consider here images as bandlimited functions with a finite number frequencies, since this is a reasonable model for restoring digital images. Moreover this model is adapted to compute convolutions and permits to include an a priori decay of the Fourier coefficients.

Based on the above considerations, we propose the following discrete regularization functional, If u ∈ B M , then we define

J d A (u) = 1 M 2 0≤r,l<M A(D)u rN M , lN M . (4.1)
If A(iω) satisfies (3.2), then J A (u) and J d A are seminorms in B M and the only function andJ d A are norms in the finite dimensional quotient space B M /R, hence they are equivalent. Notice that if u ∈ B M and we define

u ∈ B M such that J A (u) = 0 (resp. such that J d A (u) = 0) is u = constant. Thus J A (u)
J d,k A by replacing M by k in (4.1), then J d,k A (u) → J A (u) as k → ∞.
Unless we intend to zoom and restore the images, we take M = N , where N 2 is the number of data.

Assume that the input data (measurements) consist of N 2 samples {z(ξ k )} k , and let σ, β > 0, and G ∈ ℓ ∞ (Z 2 ) be a discrete, positive, normalized convolution kernel such that G(r, l) ≥ 0 and (r,l)∈Z 2 G(r, l) = 1. Then we propose to minimize the functional min u∈BN J β,A (u) := 0≤r,l<N

β 2 + |A(D)u (r, l)| 2 , subject to G * |∆ Ξ (h * u) -z| 2 (ξ k ) ≤ σ2 ∀ξ k ∈ Ξ, and
0≤r,l<N u (r, l) = ξ k ∈Ξ w k z(ξ k ) =: z w . (4.2)
To avoid the non differentiability of J d A at 0, it is a common approximation to use the regularizer J β,A instead.

Therefore we will minimize J β,A (u) on B N subject to the family of constraints in (4.2). Notice that we have incorporated the image acquisition model (1.1) as a set of local constraints. The convolution of G and v ∈ ℓ ∞ (Ξ) is defined in the usual way by imposing an arbitrary regular grid structure on Ξ, i.e.

(G * v)(ξ k ) = l∈Z 2 G(l -k)v(ξ l ).
In the case of perturbed sampling this regular grid structure may be determined by the original unperturbed grid, otherwise it may be based on a nearest neighbor computation.

Notice that we have used the value σ > 0 as an estimate of the standard deviation of the noise. We will make tests with σ = σ and also with values of σ different from σ. Also, the effective support of G must permit the statistical estimation of the variance of the noise. In Section 7 we will come back to the noise estimation issue and the choice of σ.

In the rest of the paper, we assume that the blurring kernel h satisfies

h ∈ L 2 (Ω N ), supp ĥ ⊆ [- M 2 , M 2
] 2 , and ĥ(0, 0) = 1. (

If u ∈ B M , then we can compute h * u using the Fourier transform h * u (p, q) = ĥ (p, q) û (p, q). The last equality constraint in (4.2) fixes the global mean of u to be the mean of the samples z, weighted by the areas w k (see (2.6)). This constraint is necessary to assure the uniqueness of the solution, since the data fitting is provided only by inequality constraints, and the solution may be undetermined up to a constant (in the kernel of the regularizer) in areas where the variance of z is smaller than σ 2 [START_REF] Bertalmío | TV based image restoration with local constraints[END_REF]. The details of the uniqueness proof are given in Section 5. Now, our purpose is to prove that the constrained formulation of (4.2) can be solved using Uzawa's method once we guarantee that the assumptions of Uzawa's method [START_REF] Faurre | Analyse numérique. Notes d'optimisation[END_REF] hold. But before that, we comment on a improved discretization for approximating J A (u).

An improved discretization formula.

In this Section we follow the proposal made by Moisan in [START_REF] Moisan | How to discretize the Total Variation of an image?[END_REF] to improve the discretization of the total variation formula. The basic observation is that, the computation of the Dirichlet integral in [0, N ) 2 ( [0,N ) 2 |∇u| 2 ) cannot be done accurately unless we previously zoom the image u by a factor of two.

The same argument applies to the case of ΩN |A(D)u| 2 , but not for ΩN β 2 + |A(D)u| 2 with β ≥ 0. Indeed, in the last case, an exact computation would involve an infinite number of samples. However, as Moisan has shown for the TV case [START_REF] Moisan | How to discretize the Total Variation of an image?[END_REF], doubling the number of variables leads to a good approximation of the above integral, being a good compromise between precision and algorithmic efficiency.

Let us analyze the implications of Moisan's discretization for our restoration model. To do so, we need to introduce some notation. For each M ∈ N, we denote by X M the Euclidean space R M×M . The Euclidean scalar product and the norm in X M will be denoted by •, • XM and • XM , respectively, but in absence of ambiguities we will omit the subindex. X M represents the space of images B M sampled in the regular grid {0, . . . , M -1} × {0, . . . , M -1} (or given by its Fourier coefficients û ∈ ℓ 2 (I M ). Let us introduce the operator

P : X N → X 2N , P {u(k, l)} = ū r 2 , l 2
r,l∈{0,...,2N -1} , where {u(r, l)} r,l∈{0,...,N -1} ∈ X N and ū is the function of B N defined by the samples {u(r, l)} r,l∈{0,...,N -1} . Observe that P u(p, q) = 1 IN (p, q)û(p, q). So we may consider the operator A(D) as acting on X 2N or as acting on X N . Notice that if u ∈ X N we may write A(D)P u = P A(D)u. From now on we will avoid (except in ambiguous cases) the use of subindexes to specify the function spaces of norms and scalar products, and the function space should be clear from the context. Thus, our final restoration model is

min u∈XN K β,A (u) := 0≤r,l<2N β 2 + |A(D)P u (r, l)| 2 , subject to G * |∆ Ξ (h * u) -z| 2 (ξ k ) ≤ σ2 ∀ξ k ∈ Ξ , and
0≤r,l<N u (r, l) = ξ k ∈Ξ w k z(ξ k ) =: z w . (4.4)
5. The well-posedness of the model and its numerical solution. Proposition 5.1. Assume that (4.3) holds. Then there exists a unique minimum u ∈ X N of (4.4).

Proof. Let u m be a minimizing sequence of (4.4). Since A(D)P u m is bounded in X 2N , and ω = 0 is the only vanishing frequency for A(iω) we deduce that v m := u mu m (0, 0) is bounded in X N . Now, since u m (0, 0) is constrained to be z ω , we have that u m is bounded in X N . By extracting a subsequence, if necessary, we may assume that u n → u. It is immediate to see that u satisfies the constraints. Since K β,A is lower semicontinuous, we have that u is a minimum of (4.4). Now, let u 1 , u 2 be two minima of (4.4). If A(D)P u 1 = A(D)P u 2 , letting u = u1+u2 2

, then the strict convexity of K β,A proves that K β,A (u) < inf u∈XN K β,A , a contradiction. Thus A(D)P u 1 = A(D)P u 2 and we have uniqueness modulo constants, i.e., u 1u 2 = c for some c ∈ R. Since û1 (0, 0) = û2 (0, 0) we deduce that c = 0, and therefore, u 1 = u 2 .

Remark 4. Proposition 5.1 is also true if instead of assuming the average constraint in (4.4) we assume that inf c∈R G * (zc) 2 > σ2 . This can be proved as in [START_REF] Bertalmío | TV based image restoration with local constraints[END_REF][START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF]. But in that case, we should also use a different (gradient descent based) algorithm to minimize (4.4) as described in [START_REF] Bertalmío | TV based image restoration with local constraints[END_REF].

From now on, we assume that the constraints are qualified, that is there is u

∈ X N such that û(0, 0) = z w and G * |∆ Ξ (h * u) -z| 2 (ξ k ) < σ2 , ∀ξ k ∈ Ξ, (5.1)
which implies that the set of functions satisfying the constraints is non-empty. We prove that the solution of (4.4) can be computed by adapting Uzawa's algorithm. Let µ > 0 and λ = (λ k ) N 2 k=1 ≥ 0. Define the Lagrangian function

L µ (u, {λ}) = K β,A (u) + µ (û(0, 0) -z w ) 2 + ξ k ∈Ξ λ k 2 ([G * |∆ Ξ (h * u) -z| 2 ](ξ k ) -σ2 ).
To adapt Uzawa's algorithm we need the following result which can be proved as in the proof of Proposition 5.1. Theorem 5.2. For each λ = (λ k ) N 2 k=1 ≥ 0, there is a unique solution u of

min u∈XN L µ (u, {λ}).
Proof. Since L µ (u, {λ}) is lower semicontinuous in u, it suffices to prove that any minimizing sequence u n is bounded. Since L µ (u m , {λ}) is bounded, we know that A(D)(u m ) and ( u m (0, 0)z w ) 2 are bounded. The boundedness of A(D)(u m ) implies that u mu m (0, 0) is bounded. Since ( u m (0, 0)z w ) 2 is bounded, then u m (0, 0) is also bounded.

We solve (4.4) with Uzawa's algorithm.

Algorithm II: Restoration with local constraints 1. Choose any set of values λ 0 k ≥ 0, k = 1, . . . , N 2 , and µ 0 > 0. Iterate from p = 0 until convergence of λ p the following steps: 2. With the values of λ p , µ p solve:

u p = arg min u L µ p (u, {λ p }).
(5.2)

3. Update µ and λ in the following way:

µ p+1 = µ p + 1, λ p+1 k = max(λ p k + ρ p ([G * |∆ Ξ (h * u p ) -z| 2 ](ξ k ) -σ2 ), 0) ∀ξ k , (5.3)
where 0 < ρ * ≤ ρ p ≤ ρ * . Proposition 5.3. Assume that there exists u ∈ X N such that û(0, 0) = z w and z(ξ k ) = h * u(ξ k ) ∀ξ k ∈ Ξ. Then Uzawa's algorithm converges to the solution of (4.4).

To prove Proposition 5.3 we need to reformulate problem (4.4) as

min u∈XN max λ≥0,α+,α-≥0 L(u, {λ}, α + , α -), (5.4
)

where λ = (λ k ) N 2 k=1 , α + , α -≥ 0, L(u, {λ}, α + , α -) = K β,A (u) + ξ k ∈Ξ λ k ([G * |∆ Ξ (h * u) -z| 2 ](ξ k ) -σ2 ) + α + ϕ + (u) + α -ϕ -(u),
and ϕ + (u) := û(0, 0)z w and ϕ -(u) = -û(0, 0) + z w , so that the equality constraint is written as the two inequalities ϕ + (u) ≤ 0, ϕ -(u) ≤ 0. Since, by Proposition 5.1, problem (4.4) has a solution, the classical existence result of saddle points (see [START_REF] Faurre | Analyse numérique. Notes d'optimisation[END_REF]) proves the existence of a solution of (5.4). Indeed the following result is classical and can be found, for instance, in [START_REF] Faurre | Analyse numérique. Notes d'optimisation[END_REF] (Theorems 4 and 6, pp. 59-61) or [START_REF] Ciarlet | Introduction à l'analyse numérique matricielle et à l'optimisation[END_REF] (Theorem 9.3.2).

Theorem 5.4. Assume that (5.1) holds. Let u be the solution of (4.4). Then there are ({λ}, α + , α -) ≥ 0 such that (u, {λ}, α + , α -) is a solution of (5.4), i. e., a saddle point of L(•, •, •, •). If (u, {λ}, α + , α -) is a solution of (5.4), then u is a solution of (4.4).

Since we will need it below, let us compute the gradient of K β,A (u). For any v ∈ X N we have

∇K β,A (u), v XN = A(D)P u β 2 + |A(D)P u| 2 , A(D)P v X2N = P * A(D) * A(D)P u β 2 + |A(D)P u| 2 , v

XN

for each v ∈ X N vanishing on the boundary of {0, . . . , N -1} 2 . Thus, we may write

∇K β,A (u) = A(D) * P * A(D)P u β 2 + |A(D)P u| 2 ∈ X N .
Now, we notice that P * f is just the restriction operator (subsampling operator) that considers only the samples of f ∈ X 2N in the grid {0, . . . , N -1} 2 . Finally, using this and the last two formulas, we deduce that

A(D)P u β 2 + |A(D)P u| 2 • ν {0,...,N -1} 2 = 0
where ν {0,...,N -1} 2 is the discrete normal.

As usual, we denote by v q = N 2 i,j=1 |v(i, j)| q 1/q for any v ∈ X N , 1 ≤ q < ∞. We denote v ∞ = max (i,j)∈{1,...,N } |v(i, j)|. And for simplicity in the cases where there is no ambiguity, we shall omit the subindexes for the L 2 -norm, then u 2 = u .

Proof of Proposition 5.3. Let us write Q(u) = ( u(0, 0) -z w ) 2 , R(u) = G * |∆ Ξ (h * u) -z| 2 .
To adapt the convergence proof of Uzawa's method to our case, we need to prove that (a) If U is a bounded subset of X N then there is a constant α > 0 such that

∇K β,A (u) -∇K β,A (v), u -v + µ ∇Q(u) -∇Q(v), u -v ≥ α u -v 2 for all u, v ∈ U . (b) R(u)
is Lipschitz on bounded sets of X N and (c) the sequence u p constructed in Step 2 of the above algorithm is bounded in X N . To prove (a) we use the inequality [START_REF] Trudinger | On the analyticity of generalized minimal surfaces[END_REF]

ξ β 2 + |ξ| 2 - ξ ′ β 2 + |ξ ′ | 2 , ξ -ξ ′ ≥ β 2 |ξ -ξ ′ | 2 (β 2 + |ξ| 2 + |ξ ′ | 2 ) 3/2 ∀ξ, ξ ′ ∈ R k
with k = 2 and we compute

∇K β,A (u) -∇K β,A (v), u -v ≥ (r,l) |A(D)P u(r, l)) -A(D)P v(r, l)| 2 (β 2 + |A(D)P u(r, l)| 2 + |A(D)P v(r, l)| 2 ) 3/2 ≥ α (r,l) |A(D)P u(r, l) -A(D)P v(r, l)| 2
where (r, l) ∈ {0, 1, . . . , 2N -1}, the constant α > 0 depends on the bound for U , and

∇Q(u) -∇Q(v), u -v = 2(û(0, 0) -v(0, 0)) 2 .
Then (a) follows as a consequence of the two previous inequalities.

(b) Assume that U ⊆ X N is a bounded set. Let u, u ∈ U . Since G 1 ≤ 1, we have R(u) -R(v) ≤ G 1 (h * u -z) 2 -(h * v -z) 2 ≤ 2 z ∞ h * (u -v) + h * (u + v) ∞ h * (u -v) ≤ C u -v
where C is a constant depending on the norms of h and z and on the bound for U . (c) To prove that {u p } p is bounded we observe that

L µ (u p , {λ p }) ≤ L µ (u, {λ p }), ∀u ∈ X N , (5.5)
for all p. Choosing u ∈ X N such that û(0, 0) = z w and z = ∆ Ξ (h * u), we obtain that

K β,A (u p ) + µ p Q(u p ) ≤ K β,A (u),
hence {u p } p is bounded in X N . Now, we can adapt the proof of Uzawa's method to our case (see Theorem 5 in [START_REF] Faurre | Analyse numérique. Notes d'optimisation[END_REF], Sect. 3.1). Since u p satisfies (5.5) we have ∇K β,A (u p ), u-u p +µ p ∇Q(u p ), u-u p + λ p , R(u)-R(u p ) ≥ 0 ∀u ∈ X N . (5.6)

Let u * be the solution of problem (4.4). Since, by Theorem 5.4, we have

L(u * , {λ}, α + , α -) ≤ L(u, {λ}, α + , α -) ∀u ∈ X N ,
we also have

∇K β,A (u * ), u -u * + λ, R(u) -R(u * ) +α + (ϕ + (u) -ϕ + (u * )) + α -(ϕ -(u) -ϕ -(u * )) ≥ 0 (5.7)
for all u ∈ X N . Since u * is a solution of (4.4), we have that ∇Q(u * ) = 0, and we can add µ p ∇Q(u * ), u-u * to the inequality (5.7). Taking u = u p in this form of the second inequalities, and u = u * in (5.6) and adding both of them we obtain

∇K β,A (u p ) -∇K β,A (u * ), u p -u * + µ p ∇Q(u p ) -Q(u * ), u p -u * -α + (ϕ + (u p ) -ϕ + (u * )) -α -(ϕ -(u p ) -ϕ -(u * )) + λ p -λ, R(u p ) -R(u * ) ≤ 0 Since ∇K β,A (u p ) -∇K β,A (u * ), u p -u * ≥ α A(D)P u p -A(D)P u * 2 and µ p ∇Q(u p ) -Q(u * ), u p -u * -α + (ϕ + (u p ) -ϕ + (u * )) -α -(ϕ -(u p ) -ϕ -(u * )) = 2µ p ( u p (0, 0) -u * (0, 0)) 2 -(α + -α -)( u p (0, 0) -u * (0, 0)) ≥ µ p ( u p (0, 0) -u * (0, 0)) 2
for p large enough, we have

λ p -λ, R(u p ) -R(u * ) ≤ -α A(D)P u p -A(D)P u * 2 -µ p ( u p (0, 0) -u * (0, 0)) 2 ≤ -α 0 u p -u * .
(5.8) Now, the proof follows in a standard way. Let us give the details for the sake of completeness. Using (5.3), we have

λ p+1 -λ ≤ λ p -λ + ρ p (R(u p ) -R(u * )) .
Taking squares, we have

λ p+1 -λ 2 ≤ λ p -λ 2 + 2ρ p λ p -λ, R(u p ) -R(u * ) + ρ 2 p R(u p ) -R(u) 2 .
Using (5.8) and (b), we have

λ p+1 -λ 2 ≤ λ p -λ 2 -2α 0 ρ p u p -u * 2 + ρ 2 p L 2 u p -u * 2 .
for some L > 0. If we choose ρ p such that

2α 0 ρ p -L 2 ρ 2 p ≥ γ > 0,
that is, 0 < ρ * ≤ ρ p ≤ ρ * , we have

λ p+1 -λ 2 ≤ λ p -λ 2 -γ u p -u * 2 .
Then we deduce that λ pλ is decreasing and, thus, has a limit ℓ ≥ 0. Then letting p → ∞ we have that u pu * → 0.

6. A Quasi-Newton algorithm for the solution of (5.2). The purpose of this Section is to explain the algorithm used to solve problem (5.2) in Algorithm II. For convenience, let us denote the convolution and irregular sampling operators, as Su = ∆ Ξ (h * u) for any u ∈ X N .

Observe that the Euler-Lagrange equation corresponding to (5.2) is

A(D) * P * A(D)P u β 2 + |A(D)P u| 2 + 2µ(û(0, 0) -z w ) + S * (G * λ)(Su -z) = 0. (6.1)
To shorten our expressions, let us define the following operators:

Mu = ( √ G * λ)Su, so that M * u = S * ( √ G * λu), N u = û(0, 0), b = S * ((G * λ)z) + 2µz w , A[v](u) = A(D) * P * A(D)P u β 2 + |A(D)P v| 2
and

T [v](u) = A[v](u) + 2µN u + M * Mu.
where u, v ∈ X N . We want to solve (6.1) with a fixed point iteration:

A[u t ](u t+1 ) + 2µN u t+1 + M * Mu t+1 = b. (6.
2)

The rest of this Section is devoted to show that such a fixed point algorithm converges to the minimizer of (5.2). Notice that a proof based on convex analysis (the half-quadratic regularization approach) can be found in [START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF] in the continuous case, or in [START_REF] Charbonnier | Deterministic edgepreserving regularization in computed imaging[END_REF][START_REF] Aubert | Mathematical Problems in Image Processing[END_REF] for the discrete case. Further analysis can be found in [START_REF] Acar | Analysis of total variation penalty methods for ill-posed problems[END_REF][START_REF] Dobson | Convergence of an iterative method for total variation denoising[END_REF][START_REF] Vogel | Iterative methods for total variation denoising[END_REF]. Here, we will extend the proof proposed by Chan and Mulet in [START_REF] Chan | On the convergence of the Lagged Diffusivity Fixed Point Method in total variation image restoration[END_REF]. The advantage of such an approach is that only basic algebra is needed. Moreover, the linear convergence rate of this algorithm can be shown explicitly. The difference with the approach in [START_REF] Chan | On the convergence of the Lagged Diffusivity Fixed Point Method in total variation image restoration[END_REF] relies on the fact that the operator M is not assumed to be invertible. In our case the presence of the mean constraint µN permits to prove the same convergence result as in [START_REF] Chan | On the convergence of the Lagged Diffusivity Fixed Point Method in total variation image restoration[END_REF], without an invertibility hypothesis on M.

Remark 5. Note that computing M * M in (6.2) entails the computation of an operator with a Toeplitz structure. This is efficiently computed in O(N 2 log 2 N 2 ) steps, as in the ACT algorithm (Section 2.1).

Remark 6. As mentioned in Remark 3, we may write A(D)u = DΨu. This permits to use the change of variables v = Ψu and write the constrained restoration problem (4.4) in terms of v. Then the regularizer coincides with the total variation applied to v. The solution of (4.4) can be recovered as u = Ψ -1 v. This represents no change in the developments of this paper and we keep the notation A(D)u.

6.1. Existence of u t and its boundedness. The sequence u t will be defined iteratively using (6.2).

Proposition 6.1. The equation (6.2) has a unique solution u t+1 ∈ X N which is the minimizer of

E(u) = A(D)P u (β 2 + |A(D)P u t | 2 ) 1/4 2 + 1 2 Mu -z ′ 2 + µ N u -z w 2 (6.3)
where

A(D)P u (β 2 +|A(D)P u t | 2 ) 1/4 2 = (r,l) A(D)P u(r,l) (β 2 +|A(D)P u t (r,l)| 2 ) 1/4 and z ′ = √ G * λz.
Proof. It is standard that (6.3) admits a unique solution u t+1 ∈ X N . Moreover, (6.2) is the Euler-Lagrange equation associated to (6.3) and solutions of (6.2) are minimizers of (6.3). Proposition 6.2. (i) There exists K 0 > 0 such that A(D)P u t ≤ K 0 .

(ii) T [u t ] is a bounded coercive operator. Indeed we have

T [u t ]u, u ≥ α u 2 (6.4)
for some α > 0 independent of t.

(iii) The sequence u t is uniformly bounded. Proof. Since u t+1 is a minimizer of E(u), we have E(u t+1 ) ≤ E(0) = 1 2 z ′ 2 + µ z w 2 , and thus:

1 (β 2 + |A(D)P u t | 2 ) 1/4 A(D)P u t+1 2 ≤ E(0) = 1 2 z ′ 2 + µ z w 2 .
We have

A(D)P u t+1 (β 2 + |A(D)P u t | 2 ) 1/4 2 ≥ A(D)P u t+1 2 β 2 + |A(D)P u t | 2 ∞ Thus A(D)P u t+1 2 ≤ E(0) β 2 + |A(D)P u t | 2 ∞ ≤ E(0) β 2 + A(D)P u t 2
∞ . But since we deal with finite dimensional spaces, there exists L > 0 which does not depend on u t such that A(D)P u t ∞ ≤ L A(D)P u t . Hence we deduce that

A(D)P u t+1 2 ≤ E(0) β 2 + L 2 A(D)P u t 2
Assume that A(D)P u t ≤ K. Using (6.1), to get that A(D)P u t+1 ≤ K, it is sufficient to choose K > 0 large enough so that

E(0) β 2 + L 2 K 2 ≤ K 2 .
(ii) The boundedness of T [u t ] is immediate and we omit its proof. Let us prove that T t is a coercive operator. Using the bounds in Step (i), we have 1

β 2 + |A(D)P u t | 2
A(D)P u, A(D)P u ≥ 1

β 2 + L 2 K 2 0 A(D)P u 2 .
Hence

A[u t ]u, u ≥ 1 β 2 + L 2 K 2 0 A(D)P u 2 ≥ α 0 u 0 2 ,
where we wrote u = u 0 + c, with c = N u and N u 0 = 0, and we used the fact that A(D)P is a linear operator whose kernel are the constants. Using A[u t ]c = 0 and u 0 , c = 0, we get

T [u t ](u 0 + c), u 0 + c ≥ A[u t ](u 0 + c), u 0 + c + 2µ N (u 0 + c), u 0 + c ≥ A[u t ]u 0 , u 0 + 2µc 2 ≥ α 0 u 0 2 + 2µc 2 .
Thus, we deduce (6.4).

(iii) From (6.4), we know that T [u t ]u t+1 , u t+1 ≥ α u t+1 2 . But from (6.1), we know that T [u t ]u t+1 , u t+1 = b, u t+1 ≤ b u t+1 . We deduce that u t+1 ≤ b α .

6.2. Convergence of the fixed point algorithm. For simplicity, given λ = (λ k ) N 2 k=1 , we write L(u) = L µ (u, {λ}). Recall that

∇ u L(u) = T [u](u) -2µz w -M * z ′ .
Let us finally define

G(v, u) = L(u) + v -u, ∇ u L(u) + 1 2 v -u, T [u](v -u)
Proposition 6.3. The following inequality holds for any u, v ∈ X N :

L(v) ≤ G(v, u) (6.5)
Proof. We follow the proof in [START_REF] Chan | On the convergence of the Lagged Diffusivity Fixed Point Method in total variation image restoration[END_REF]. Since

G(v, u) -L(v) = L(u) -L(v) + v -u, ∇ u L(u) + 1 2 v -u, T [u](v -u) ,
standard computations lead to

G(v, u) -L(v) = (r,l) a -a + 1 2a (a 2 -a 2 ) with a = β 2 + |A(D)P u(r, l)| 2 and a = β 2 + |A(D)P v(r, l)| 2 where (r, l) ∈ {0, 1, . . . , 2N -1}. Since a, a > 0, and a -a + 1 2a (a 2 -a 2 ) = (a-a) 2 2a ≥ 0, we have that G(v, u) -L(v) ≥ 0. Proposition 6.4. (i)
The function u t+1 defined by (6.2) is such that:

u t+1 = argmin v G(v, u t ) (6.6)
i.e.: 0 = ∇ u L(u t ) + T [u t ](u t+1u t ).

(ii) We have lim t→+∞ u t+1u t = 0.

Proof. (i) Let us denote by ū = argmin v G(v, u t ). We thus have 0 = ∇ u L(u t ) + T [u t ](ūu t ). And this last equation is precisely equation (6.2), which implies that ū = u t+1 .

(ii) From (6.5) and (6.6), we have L(u t+1 ) ≤ G(u t+1 , u t ) ≤ G(u t , u t ) ≤ L(u t ), i.e. (L(u t )) is decreasing. Now, from (6.5) and (6.6), we have:

L(u t+1 ) ≤ G(u t+1 , u t ) = L(u t ) + u t+1 -u t , ∇ u L(u t ) + 1 2 u t+1 -u t , T [u t ](u t+1 -u t ) = L(u t ) - 1 2 u t+1 -u t , T [u t ](u t+1 -u t )
Using (6.4), we deduce:

1 2 α u t+1 -u t 2 ≤ 1 2 u t+1 -u t , T [u t ](u t+1 -u t ) ≤ L(u t ) - L(u t+1
) and (ii) follows.

We are now in position to state a convergence result. Theorem 6.5. The sequence u t defined by (6.2) converges to the solution of (5.2).

not desirable to compute the solution that (absolutely) minimizes the TV but to keep a solution with a slightly higher TV value in order to avoid the loss of textures.

As a consequence, to avoid this degradation, the rule of thumb has been ever since to remove less noise than noise is actually present in the image. Gilboa studied this in [START_REF] Aujol | Constrained and SNR-based Solutions for TV-Hilbert Space Image Denoising[END_REF] and concluded that in terms of SNR the optimal selection of σ is between 0.7 and 0.8 times the value of σ. In what follows we will modify the constraints to account for this change

G * |∆ Ξ (h * u) -z| 2 (ξ k ) ≤ σ2 ∀ξ k ∈ Ξ ,
where σ < σ. That is, in order to keep more texture in u we do not remove all noise and we could write z = S(u 0 + S -1 (nn)) + n where n is a noise with variance σ2 and we identify u 0 + S -1 (nn) with u.

Motivated by this observation and by the fact that, due to noise, there is always some uncertainty in the surroundings of ∂V , we will avoid the computational overhead of getting exactly to ∂V by stopping the algorithm as soon as the solution is close to it. This is the notion behind the band constraint

(1 -α)σ 2 ≤ G * |∆ Ξ (h * u) -z| 2 (ξ k ) ≤ (1 + α)σ 2 ∀ξ k ∈ Ξ, ( 7.1)
with α > 0. The constraint described by equation (7.1) is clearly non-convex, and therefore it cannot be integrated in the method presented here. But since Uzawa's algorithm always pushes the solution towards the boundary of the feasible set, then (at least in practice) it can be used to stop Uzawa's loop by testing if (7.1) is fulfilled.

In what follows we will see that even considering relaxed constraints like (7.1), imposing all of the local constraints simultaneously is not statistically correct, since each constraint's estimator behaves as a random variable. Then we will see how this relaxation of constraints is used to early stop Uzawa's iterations and how this helps to improve the efficiency of our implementation.

In our experiments, we have chosen σ = 0.8σ and α such that 0.8(1 + α) < 1.

7.1. Expected number of satisfied band constraints. Let us summarize the arguments of [START_REF] Almansa | A TV based restoration model with local constraints[END_REF] and adapt them to the case of band constraints. Each local constraint relies on a local estimate of the residual variance of the form

S G (ξ k ) = G * |∆ Ξ (h * u) -z| 2 (ξ k ) = G * |n| 2 (ξ k ) , (7.2)
where G is a Gaussian or uniform window centered at the interest point and n k denotes a zero mean Gaussian noise and variance σ2 (recall that we are going to remove only a noise of variance σ2 < σ 2 ). Since S G is a random variable itself, the number of satisfied constraints is estimated by the probability

P [S G ≤ (1 + α)σ 2 ] in the case of the ball shaped constraint, or P [(1 -α)σ 2 ≤ S G ≤ (1 + α)σ 2 ] for the band constraint.
Observe that a constraint of the type S G = σ2 (that is in practice imposed when solving exactly (4. Then, for the case of the ball constraint, the expected number of satisfied constraints is the number of pixels times the following probability

P [S I ≤ (1 + α)σ 2 ] = P 1 |I| k∈I n2 k ≤ (1 + α)σ 2 = P k∈I nk σ 2 ≤ (1 + α)|I| .
Notice that in the rightmost equation k∈I nk σ 2 is a sum of |I| squared normalized Gaussian random variables, so it follows a chi-square distribution with |I| degrees of freedom (χ 2 (|I|)). And the probability can be computed using the incomplete gamma function

Γ(a, x) = ∞ x t a-1 e -t dt P [S I ≤ (1 + α)σ 2 ] = P χ 2 (|I|) ≤ (1 + α)|I| = Γ (1 + α)|I| 2 , |I| 2 .
In the case of the band constraint the expected number of satisfied constraints N (α, r) is the number of pixels times the following probability

P (1 -α)σ 2 ≤ S I ≤ (1 + α)σ 2 = Γ (1 + α)|I| 2 , |I| 2 -Γ (1 -α)|I| 2 , |I| 2 . (7.3)
Equation ( 7.3) expresses the expected proportion of satisfied constraints as a function of the radius of the disk r (|I| = πr 2 ) and the width of the band α. We plot in Figure 7.2 this function, for different values of α to give an intuition of its behavior. Notice that the expected number of satisfied constraints decreases as we reduce the band width α or the radius r. This permits to determine one parameter as a function of the other two, i.e. using a disk of radius r = 13 (or a Gaussian with standard deviation 7.0) and defining a band of width 0.2σ 2 (α = 0.1) gives 89% of satisfied constraints. In practice, either we specify α, r and then the expected number of satisfied constraints is N (α, r), or we give α and N α and we compute r so that N (α, r) = N α . We have taken the second option in the experiments displayed in Section 8. Remark 7. Observe that in (7.2) the estimation of the noise variance corresponds to the case when the mean of the random variable is zero. Indeed, the global mean is enforced to be zero in (4.4). We should also impose that the local means are zero with a new set of constraints, otherwise S G will be an overestimation of the noise variance. Adding the local mean constraint (4.4) adapts to the formalism developed in this paper. But to avoid the computational overhead of its implementation, and since the overestimation plays in favor of the relaxation arguments presented earlier in this Section, we will not include it in the present formulation. The first sequence shows the error evolution along the Quasi-Newton iterations, without applying the stopping condition for the CG algorithm. The second sequence was obtained stopping the CG with the empirical bound, observe that the two sequences are indistinguishable. The third sequence depicts the estimated error used to effectively stop the Quasi-Newton iterations, as soon as the desired error is achieved.

ξ k ∈Ξ |G * (∆ Ξ (h * u) -z)| 2 (ξ k ) = 0 in
Table 7.1 Relation between the width of the constraint band and the restoration time. Reducing the width of the band increases the computational cost of the restoration algorithm. In the table the value of η was selected according to the rule η = 7.2. Practical stopping conditions for an efficient implementation. Using (7.3) we may derive a practical rule to stop Uzawa's loop. Indeed, the user specifies α and the proportion of constraints N α /N that must lie within the band of width α, and the algorithm deduces the radius r of the kernel G such that N (α, r) = N α . Then we iterate the Uzawa's loop in Algorithm II until the number of pixels that satisfy the constraint (7.3) is at least N α .

The truncation error of the Quasi-Newton has three sources: (i) truncation of the Quasi-Newton iterations themselves, (ii) truncation of the nested CG loop, and (iii) propagation of the CG error along QN iterations. Here we summarize how to estimate and control the combination of the three errors for a given (global) target error bound on the QN result u pu t p 2 ≤ η/ S . Using standard error propagation analysis [START_REF] Dahlquist | Numerical Methods[END_REF] and the knowledge that QN is at least linearly convergent we can estimate the global error determining (ii) and bound the inverses of the operators T [u t ] and their dependence on u t . In our case these bounds are estimated empirically, and the CG error is approximated by its residual. Figure 7.3 shows that this procedure is quite effective in practice. First our CG stopping condition makes the truncated QN sequence indistinguishable from the non-truncated QN sequence (i.e. the one with CG iterated until exact convergence is reached, and thus not affected by CG truncation errors (ii) and their propagation (iii)). This is because the actual QN truncation error (with either CG truncated or not) is considerably over-estimated by our error propagation analysis as shown in Figure 7.3. This also means that the desired error is achieved much faster than predicted by the our error bounds.

Let us now take into account the truncation error of the Quasi-Newton method in the determination of the band constraints. Assume that we are computing Quasi- Newton solution u p and we have controlled the errors Su p -Su t p ≤ S u p -u t p ≤ η, where u t p denotes the solution obtained at the t-th iteration of the method. Then we may erode the band by η as seen in Figure 7.1. In this way we ensure that if we stop the Quasi-Newton solution with the criterion S u pu t p ≤ η, then the truncated solution u p satisfies (7.3). Clearly η must satisfy the inequality η < 1 2 ( √ 1 + ασ -√ 1ασ) and we have taken η := 1 3 ( √ 1 + ασ -√ 1ασ), which equalizes the widths of the three bands in Figure 7.1. Figure 7.1 illustrates the band and its reduced version for a single constraint.

Finally, we notice that the computational complexity of the algorithm increases as we reduce the width of the band. Indeed, taking α → 0 makes it harder to satisfy the constraints. And to illustrate this we display in Table 7.1 the computation times corresponding to different values of α.

8. Experiments. We will test the proposed algorithm in three contexts, first in the irregular to regular sampling and denoising task, then in Subsection 8.2 we add the deconvolution, and finally in Subsection 8.3 we consider the full restoration problem, with deconvolution, denoising, and zooming.

We compare the results obtained by the ACT Algorithm [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF] (when applicable), the algorithm described in [START_REF] Almansa | Restoration and zoom of irregularly sampled, blurred and noisy images by accurate total variation minimization with local constraints[END_REF] (ACT+TV), and the proposed algorithm discussing different choices of A(D) in each context. In addition we will also consider a reformulation of our local constraint algorithm based on the minimization of the L 2norm of the pseudo-differential operator A(D). The L 2 -norm formulation is faster and leads to a linear system (which is solved directly with Conjugate Gradient, avoiding the need of a fixed point loop), and it is derived by replacing the regularizer K A,β in ( All the experiments were performed with simulated images. The perturbations ε(x) were generated according to the model (2.2), with an amplitude A = 0.88 (stan-Table 8.1 Comparison of the algorithms in the irregular to regular sampling and denosing task. These results correspond to the restoration of the image shown in Figure 8.1, corrupted by a white Gaussian noise with standard deviation σ = 1. The error column is obtained by comparing the restored image u with the ground truth u 0 , where RMSE(u, u 0 ) = ( 1N u -u 0 2 ) 1/2 . Its values evidence that all the algorithms achieve errors similar to the noise variance. In all the experiments the power of the removed noise was set to be 1 N Su -z 2 ∼ 0.908, a little below the noise level σ 2 = 1, but this allows to improve both the result's RMSE and the visual quality of the restored images.

Algorithm

Regularizer RMSE ACT [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF] dard deviations of ε(x)), and where supp ε ⊆ [-0.5 Tε , 0.5 Tε ] 2 for T ε = 10. The perturbed samples z were computed very accurately (usually 10 -8 ) by approximating the irregular sampling formula (2.3) with the transposed NFFT [START_REF] Potts | Fast Fourier transforms for nonequispaced data: A tutorial[END_REF]. Finally, we added a white noise of standard deviation σ to the irregular samples. In the experiment displayed here, we have taken σ = 1 gray levels (i.e.the noise power is 890 times smaller than the image power, SN R = 29.5dB). The perturbed image shown in Figure 8.1, which corresponds to the denoising experiment, was simulated according to this procedure.

To quantify the errors we adopt the classical root mean squared error measure RMSE(u, u 0 ) = ( 1N uu 0 2 ) 1/2 against the ground truth image (denoted as u 0 ), and the method noise for a qualitative analysis. The method noise was originally aimed at comparing denoising algorithms. It consists in subtracting the restored image to the noisy one, and studying the remaining noise. In our context assuming an image formation model like (1.1) and denoting u the image obtained by a restoration algorithm, the method noise becomes (z -Su), where z are the noisy samples and S stands for the irregular sampling and convolution operators. Since the restoration is expected to recover the original image u ≃ u 0 , the method noise should be as similar to a white noise as possible. In addition, since we would like the original image u 0 not to be altered by the restoration method, the method noise should not exceed the actual noise variance, which justifies our local constraints approach.

8.1. Denoising. Observe in Table 8.1 that the proposed algorithm outperforms (in terms of RMSE) the ACT and ACT+TV [START_REF] Almansa | A TV based restoration model with local constraints[END_REF] algorithms, in the denoising experiment. Also notice in Figure 8.2 how the method noise of the proposed algorithm retains less structure, meaning that the method removes just the noise with less alteration of the texture.

The variants of the algorithm based on the minimization of an L 2 -norm (8.1), give results that are comparable with the ones obtained with the FAR regularizer, and in both reported cases the L 2 results are indistinguishable from the ones obtained by FAR. This observation advocates for the local formulation of the constraints and it can be justified by the density of the samples in this denosing experiment, and by the lack of both deconvolution and zoom. We can conclude that, for the denoising (without either deconvolution nor zoom) the advantage of a nonlinear regularizer over the L 2 one is marginal, mostly because no spectral extrapolation is needed. However minimizing the L 2 -norm is not expected to perform well in tasks that entail a spectral extrapolation like deconvolving or zooming.

The imposition of the spectral profile produces a consistent improvement in all the cases (| 2π N ω| vs. | 2π N ω| 1.6 vs. | 2π N ω| 1.9 ), and imposing the profile corresponding to the coefficient decay of the reference image (Figure 8.1) produces the best results. We can study the performance of the spectral profiles by analyzing the frequency distribution of the errors of the restored images with respect to the reference image (GT), shown in Figure 8.3. There we can see that using the linear profile (total variation), the low frequencies are heavily penalized and most of the errors come from them, but imposing the profile corresponding to this image (| 2π N ω| 1.9 ) we reduce the errors in the low frequency range.

Finally, let us spend a word to note the impact of the practical stopping conditions proposed in Section 7, since they allow a significant speed up of the algorithm reducing the execution time from 200 sec (stopping after 50 Uzawa's iterations) to 30 sec for images of size 149 × 149 pixels. 8.2. Denoising and deconvolution. We consider in this Section the denoising and deconvolution of irregularly sampled images. For that we include in our image formation model the MTF corresponding to SPOT 5 HRG (High Resolution Geometric) satellite with Hipermode sampling [START_REF] Hagolle | How to double the spatial resolution of a pushbroom instrument[END_REF]. Shortly, Hipermode is a push-broom acquisition mode that uses two shifted bars of sensors to sample on a double-density 8.2 have a higher RMSE. deconvolution and/or zoom, due to the inability of this type of regularizers to extrapolate the spectrum. This difference is also illustrated in Figure 8.6 and in Figure 8.8 for the case of deconvolution and zoom.

Extension to zooming.

Zooming requires to interpolate and restore the image while preserving and enhancing the shapes, this can be seen as a spectrum extrapolation problem. The basic idea is to fit in as much as this is possible, the low frequency components of the restored and zoomed image to the original data, and to extrapolate the spectrum to the rest of the frequency domain by means of the regularization functional. The regularization allows to recover some high frequencies, which is indeed much more convenient than just filling them with zeros, a technique which is known to produce ringing.

Since the FAR regularizer allows to control the spectral behavior of the solution, in particular the extrapolated part, it will allow to improve the zoom results. Let us first extend the formulation (4.4) to consider the restoration of irregularly sampled

β 2 + |A(D)P u (r, l)| 2 , subject to G * |∆ Ξ (h * p * u) -z| 2 (ξ k ) ≤ σ 2 ∀ξ k ∈ Ξ and 0≤r,l<nN u (r, l) = ξ k ∈Ξ w k z(ξ k ). (8.2)
The zoomed and restored image u is a vector of size nN × nN (we recall that the size of z is N × N ), and p is a spectral projector (e.g. p = χ R or a prolate function) on a low-band region R which depends both on the MTF and the sampling set. In the context we are considering here, p will be different from zero in the frequency band corresponding to the resolution of the data [-1/2, 1/2] 2 and the constraint is saying that the data is explained by the lower frequency part h * p * u of h * u. The regularization functional K n β,A (u) penalizes the oscillations that may appear when we extrapolate the high frequencies in the spectral region [-n/2, n/2] 2 \ [-1/2, 1/2] 2 . Let us mention that, as discussed in [START_REF] Almansa | Echantillonnage, Interpolation et Détection[END_REF] in the context of regular sampling, the right choice of the spectral region R permits to reduce the aliasing effects, but we shall not consider this problem here. For us, if we want to restore and zoom the image u by a factor n, p will be different from zero on the region R = [-1/2, 1/2] 2 and zero on [-n/2, n/2] 2 \R . This is a way to impose that the restored image fits the data z at low frequencies and the high ones are extrapolated via the minimization of K n β,A (u). This minimization problem (8.2) with p(ω) = 1 IN (ω) is a direct extension of the oversampling and denoising method introduced by Malgouyres and Guichard [START_REF] Malgouyres | Edge direction preserving image zooming: a mathematical and numerical analysis[END_REF] to the more general case of irregular to regular sampling, deconvolution, denoising and oversampling.

The experiments shown here correspond to a 2X zoom, the images were simulated using the procedure described for the deconvolution case, with a filter ĥ extended with zeros up to a double size, and where the irregular sampling is performed at a double spacing. The resulting image is shown in Figure 8.9. Since the restored image and the reference image u 0 have the same size, they can be directly compared. In Figure 8.9 are shown the distorted and the reference image as well as two restorations.

Let us first comment on the stair-casing effect that is noticeable in the bottom left image of Figure 8.9. It is a common observation that the total variation introduces a stair-casing effect in the restored images, but let us point out that in our case where the derivatives are computed analytically this effect may not appear. It appears in the bottom left image of Figure 8.9 due to the poor discretization of the total variation. Notice that the same image processed with a finer approximation as proposed in Section 4.1 does not exhibit this artifact, see bottom right image in Figure 8.9. As we mentioned above, this effect is product of the coarse discretization of the total variation and was negligible in the cases of restoration without zoom. But when zooming it is important to avoid this effect since it produces un-natural looking images in spite of the fact that the RMSE errors of both images are similar. Since the spectrum is extrapolated the quality of the zoomed image is not affected by the penalty in the frequencies imposed by the regularization term. In Table 8.3 we display the results obtained with different penalization profiles. Notice that the results are better in the case of L 1 -norm than L 2 -norm.

In contrast with the previous applications, in the present case, not removing all the noise leads to some artifacts. The noise is defined over the original grid (be either regular or irregular), but any residual of the original noise becomes a low frequency colored noise in the zoomed image that is visible as artifacts. 9. Conclusions. We have proposed a model for the restoration of band limited images that considers irregular (perturbed) sampling, denoising, and deconvolution. In addition a theoretically sound discretization of this model, and an algorithm that solves the resulting discrete optimization problem efficiently using automatically optimized thresholds and stopping conditions for a prescribed output precision level, are presented.

An experimental evaluation (which actually requires all three modeling, discretization, and algorithmical elements) shows that our restoration model, and its accurate discrete approximation, actually improves the performance of previous methods in terms of both signal to noise ratio (Tables 8.1, 8.2 and 8.3) and adequacy of the statistical properties of the computed solution with respect to the assumed image acquisition model. As we can see in Figures 8.2 and 8.5 our approach permits to recover part of the structure lost with a global constraint leading to a method noise that resembles more the postulated gaussian white noise. As observed by Nikolova in [START_REF] Nikolova | Model distortions in bayesian map reconstruction[END_REF], most denoising methods tend to produce a method noise with a distribution which largely differs from the postulated gaussian white noise hypothesis.

The key ingredients that make our restoration model (Section 3) produce good results are the combination of 1. the modeling of the image acquisition system through local constraints as in [START_REF] Almansa | A TV based restoration model with local constraints[END_REF][START_REF] Almansa | Restoration and zoom of irregularly sampled, blurred and noisy images by accurate total variation minimization with local constraints[END_REF], 2. the use of a special kind of frequency adapted regularization in the spirit of Fig. 8.9. Deconvolution denoising and zoom. We display the QN+FAR results using a spectral decay like: ( 2πω N ) 1.6 , which coincides with the decay of the reference image. The stair-casing effect in (c) (RMSE: 6.003) is due to the poor discretization of the regularizer in (4.2). In (d) (RMSE: 5.957) we observe the result obtained with a finer discretization of the regularizer (4.4).

the ACT algorithm proposed by H.G. Feichtinger, K. Gröchenig, M. Rauth and T. Strohmer [START_REF] Gröchenig | Numerical and theoretical aspects of non-uniform sampling of band-limited images[END_REF], and 3. the minimization of an L 1 (instead of L 2 -norm) of this frequency-adapted pseudo-differential operator. However the L 2 -norm produces competitive results in the case of denoising without deconvolution nor zoom (see Table 8.1). All three aspects together collaborate to provide a better model of both, image acquisition and a priori image regularity, that leads to better accordance of the statistical properties of our solution with respect to the postulated models and the reference image. An accurate discrete approximation to this model (Section 4), both allows the use of convenient and efficient numerical algorithms (Section 6), and leads to wellbehaved solutions (Section 5). Finally the study of the local constraints (Section 7) leads to a trade-off between the precision of the noise estimation and the number of constraints that can be ensured, that is used to derive the stopping criterion for the different stages of the algorithm.

Experiments (Section 8) highlighting the importance of the modeling and discretization aspects described earlier are presented, namely: the quality boost provided by the local formulation of the constraints (Figure 8.2), the importance of the discretization of the regularizer (Figure 8.9) and the usefulness of the frequency adapted regularizer in favoring certain spectral behaviors of the solutions (Figure 8.2 and 8.5).

Our prospects for future work include adapting our algorithm to use splines as the underlying interpolation model, and the introduction of anti-aliasing filters to handle the case of aliased data. The exploitation of self-similarity proporties of most natural images would also be useful to increase resolution while reducing noise and aliasing effects.

Fig. 7 . 2 .

 72 Fig. 7.2. Relation between the expected number of satisfied constraints and the radius r of the disk used for the local noise estimation. The graph shows the curves computed for different widths of the band α.

Fig. 7 . 3 .

 73 Fig. 7.3. Truncation of the Quasi-Newton method. The first sequence shows the error evolution along the Quasi-Newton iterations, without applying the stopping condition for the CG algorithm.The second sequence was obtained stopping the CG with the empirical bound, observe that the two sequences are indistinguishable. The third sequence depicts the estimated error used to effectively stop the Quasi-Newton iterations, as soon as the desired error is achieved.

  ) = 0.06 s -1.9 Reference image: GT Modelled coeff. decay: f(s) (b) Mean Fourier coefficient amplitude for the reference image. (c) Noisy and irregularly sampled image.

Fig. 8 . 1 .

 81 Fig. 8.1. Image used in the denoising experiment. At left: the reference image is non-aliased, and it has 149 × 149 pixels. It was multiplied by a smooth window vanishing on the borders in order to avoid periodization artifacts (not shown). In this image the mean Fourier coefficient amplitudes decay like (2πω/N ) -1.9 (center). At right, is shown the perturbed sampled image. The perturbations have an amplitude of A = 0.88 pixels (standard deviation of ε(x)), and were simulated according to (2.2) as colored noise with the spectral content inside [-0.5/Tε, 0.5/Tε] 2 for Tε = 10.

  4.4) by L A (u) := 0≤r,l≤2N |A(D)P u(r, l)| 2 . (8.1)

  Fig. 8.2. Method noise of the different algorithms in the experiment of irregular to regular sampling plus denoising. The images in the first row display the method noise (z -Su) for different methods, less visible structure indicates a better reconstruction. The images display the grayscale range [-3, 3] scaled to [0, 255] (the full grayscale range of the image z is [0, 155]). In the second row we show the corresponding Fourier transforms, the spectrum highlights the structures that are barely visible in the spatial domain.

  norm of gradient L 2 -norm of A(D) (1.9) (b) Errors in the restored images.

Fig. 8 . 3 .

 83 Fig. 8.3. Denoising experiments with different selections of spectral profiles A(iω). Using the spectral profile that fits the model for the ground truth image, reduces the error, especially in the mid-low frequencies.

Fig. 8 . 4 .Fig. 8 . 5 .

 8485 Fig. 8.4. Image used in the deconvolution experiment. The Fourier coefficients' amplitudes in the reference image (257 × 257 pixels) decay as (2πω/N ) -1.1 .

Fig. 8 . 6 .

 86 Fig.8.[START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF]. Two details of images deconvolved using total variation and L 2 -norm of the gradient as regularizers, both with local and global constraints. From left to right is shown the result of the proposed algorithm (4.4), ACT+TV[START_REF] Almansa | Restoration and zoom of irregularly sampled, blurred and noisy images by accurate total variation minimization with local constraints[END_REF] (global constraint), the L 2 -norm with local constraints and with a global constraint. Note that the last two results are more blurred that the total variation ones, and as seen in Table8.2 have a higher RMSE.

2 (

 2 |2πω/N| 1.6 |2πω/N| 0.4 ∝ |2πω/N| 0.4 if |2πω/Ν| ≥ π/norm of gradient (b) Errors in the restored images.

Fig. 8 . 7 .

 87 Fig. 8.7. Deconvolution experiments with different selections of spectral profiles A(iω). Using the spectral profile that fits the model for the ground truth image, reduces the error.

Fig. 8 . 8 .

 88 Fig. 8.8. Deconvolution and zoom using TV and the L 2 -norm of the gradient as regularizers.The images show a detail of the zoomed images and its spectrum. Note that the L 2 regularizer recovers the low frequencies but is unable to extrapolate to higher ones.

4 K 4 L 2 -

 442 β,A (u) eq (4.4) 5.957 QN+TV |A(iω)| = | 2π N ω| 1 5.965 QN+FAR |A(iω)| ∼ | 2π N ω| 0 for |ω| > N 4 5.996 and |A(iω)| = | 2π N ω| 1 for |ω| <= N norm of the gradient L A (u) eq. (8.1) 6.270 L 2 -norm of A(D)u with |A(iω)| ∼ | 2π N ω| 1.6 for |ω| ≥ N 4 6.270 and |A(iω)| = | 2π N ω| 1 for |ω| <= N 4

  All the reported times correspond to experiments ran on a 1.6Ghz CPU restoring a 256 × 256 pixels image.

	1 3 ( √ 1 -ασ), and σ = 1. Band Width 1 + ασ -√ Effective Band Total
	parameter: α 0.60	p	(1 -α)σ + η 0.84	p (1 + α)σ -η 1.05	running time 3 min 3 sec
	0.34		0.93	1.04	3 min 37 sec
	0.22		0.96	1.03	4 min 58 sec
	0.10		0.98	1.02	18 min

Table 8 . 3

 83 Deconvolution denoising and 2X zoom experiments.

	Algorithm	Regularizer	RMSE
	ACT+TV [5]		6.236
	QN+FAR |A(iω)| ∼ | 2π N ω| 1.6 for |ω| > N

Acknowledgement. Vicent Caselles acknowledges partial support by the PNPGC project, reference MTM2006-14836. Gabriele Facciolo and Andrés Almansa acknowledge support by the ALFA project, No. AML//19.0902/97/0666/II-0366-FA. Andrés Almansa acknowledges support from CNES contracts R-S06/ OT-001-040 and R-S06/ OT-004-024 and by ECOS-Sud project U06E01. Andrés Almansa and Jean-François Aujol acknowledge support by the French "Agence Nationale de la Recherche" (ANR), under grant FREEDOM (ANR07-JCJC-0048-01), "Films, REstauration Et DOnnées Manquantes".

The relaxed optimization problem allows solutions within a band. From the depicted sequence of solutions u t p (black dots) the ones that are truly satisfying the constraints are those inside the dotted ring, since it accounts for the approximation error η. S stands for the irregular sampling and convolution operators, Su t p is also a sequence defined over the irregular samples and not over the solution space.

Proof. From Proposition 6.2.(iii), we know that u t is uniformly bounded. There exists u such that we can extract a convergent subsequence, which we still denote by u t , with u t → u as t → +∞. From Proposition 6.4.(ii) , we know that u t+1 is also convergent and u t+1 → u as t → +∞. Letting t → +∞ in equation ( 6.2), we deduce that u is the solution of (6.1) (and thus of (5.2)). By uniqueness of the solution of (5.2), we conclude that the whole sequence u t converges to u, solution of (5.2).

We end this Section by stating a result about the convergence rate. We denote by ũ the solution of Problem (5.2). We use the following notations:

and

where λ min (M ) denotes the smallest eigenvalue of the matrix M ; in particular if M is positive definite then λ min (M ) > 0 Proposition 6.6.

1. L(u t+1 ) -L(ũ) ≤ γ t (L(u t ) -L(ũ)). 2. η < 1 and 0 ≤ γ t ≤ η, for t sufficiently large. In particular, L(u t ) has a linear convergence rate of at most η. 3. u t is r-linearly convergent with a convergence rate of at most √ η. Proof. We refer the interested reader to the proof of Theorem 6.1 in [START_REF] Chan | On the convergence of the Lagged Diffusivity Fixed Point Method in total variation image restoration[END_REF] which can easily be extended to our case.

7. Band constraints and stopping conditions. Coming back to the optimization problem (4.4), since both the functional K β,A (u) and the constraints are convex, and the constraint's feasible set V does not contain the absolute minimun of K β,A (u), then the solution lies in the boundary of V . In practice, computing a solution in ∂V is not only computationally too expensive (due to the size of the problems), but also unnecessary because of the noise. Since we rely on noise estimates that have a certain accuracy, exceeding this accuracy in the data fitting is useless (as we will see later in this Section). Moreover, as it has been observed in all numerical experiments [START_REF] Blanc-Féraud | Nonlinear image processing: modeling and fast algorithm for regularization with edge detection[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF][START_REF] Chan | Continuation method for total variation denoising problems[END_REF][START_REF] Chan | On the convergence of the Lagged Diffusivity Fixed Point Method in total variation image restoration[END_REF][START_REF] Charbonnier | Deterministic edgepreserving regularization in computed imaging[END_REF][START_REF] Durand | Image Deblurring, Spectrum Interpolation and Application to Satellite Imaging[END_REF][START_REF] Geman | Constrained Image Restoration and Recovery of Discontinuities[END_REF][START_REF] Aujol | Constrained and SNR-based Solutions for TV-Hilbert Space Image Denoising[END_REF][START_REF] Gilboa | Variational denoising of partly-textured images by spatially varying constraints[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Vogel | Iterative methods for total variation denoising[END_REF], using total variation as regularizer in denoising or restoration generally carries some loss of texture and it is

where α 1 = 3.73 , β 1 = 1.75, and where sinc π (x) = sin(πx) πx if x = 0 and sinc π (x) = 1 if x = 0. This function has its first zero crossing at frequency N/2 = 1/2Hz, while at frequency 1/4Hz the power of the MTF is only 1%, meaning that outside the spectral support [-N/2, N/2] 2 (or [-1/2, 1/2] 2 Hz) of the filter there is almost no information. To simulate a blurred, noisy and irregularly sampled image with N × N samples we start by applying the MTF ĥ to the image u 0 ∈ X N . The samples z (Figure 8.4) are then obtained by sampling h * u 0 on the irregular grid and adding the noise.

In Figure 8.5 and Table 8.2 we report the results of the restoration experiments corresponding to deconvolution with the Hipermode MTF. The small variability in the results is an indicator that the MTF ĥ makes the restoration harder per-se, this is also confirmed by restoring a regularly sampled image and observing that in that case the error is similar to the irregular one (see Table 8.2). Anyway, the proposed algorithm outperforms ACT+TV [START_REF] Almansa | Restoration and zoom of irregularly sampled, blurred and noisy images by accurate total variation minimization with local constraints[END_REF] mainly thanks to the local formulation of the constraints (see Figure 8.6).

We observed that manipulating the spectral profile A(D) does not produce improvements consistent with the denoising case (see Figure 8.7). This is due, on one hand, to the fact that the reference image has a spectral decay different from the target image class (| 2π N ω| 1.6 ). Indeed the Fourier coefficients of the reference image decay as (| 2π N ω| 1.1 ), which explains the good performance of the total variation (see Table 8.2). On the other hand, in Figure 8.7(b) we see that the total variation controls the error in the low frequency range, but in the high frequency range it is too conservative and does not promote the spectral extrapolation.

The previous observation motivates the following experiment, building a profile that combines, the decay of the total variation in the low frequency range, with a decay similar to | 2π N ω| 0.4 in the high frequency range. The result of this experiment is shown in Table 8.2 and its profile is depicted in Figure 8.7(b), there we can confirm the desired effect.