Kévin Polisano
email: kevin.polisano@univ-grenoble-alpes.fr.

Laurent Condat
email: laurent.condat@gipsa-lab.grenoble-inp.fr.

Marianne Clausel

Valérie Érie Perrier

K Évin Polisano

A convex approach to superresolution and regularization of lines in images

Keywords: Super-resolution, sparse recovery, convex optimization, line detection, splitting method, spectral estimation AMS subject classifications. 94A08, 94A12, 94A20, 47N10, 90C22, 90C25

We present a new convex formulation for the problem of recovering lines in degraded images. Following the recent paradigm of super-resolution, we formulate a dedicated atomic norm penalty and we solve this optimization problem by means of a primal-dual algorithm. This parsimonious model enables the reconstruction of lines from lowpass measurements, even in presence of a large amount of noise or blur. Furthermore, a Prony method performed on rows and columns of the restored image, provides a spectral estimation of the line parameters, with subpixel accuracy.

1. Introduction. Many restoration or reconstruction imaging problems are illposed and must be regularized. So, they can be formulated as convex optimization problems formed by the combination of a data fidelity term with a norm-based regularizer. Typically, given the data y = Ax + , for some unknown image x to estimate, some known observation operator A and some noise , one aims at solving a problem like:

Find x ∈ arg min x 1 2
Axy 2 + λR(x) ,

where λ controls the tradeoff between data fidelity and regularization and R is a convex regularization functional. R can be chosen to promote some kind of smoothness.

The classical Tikhonov regularizer R(x) = ∇x 2 2 generally makes the problem easy to solve, but yields over-smoothing of the textures and edges in the recovered image x. A popular and better regularizer is the total variation R(x) = ∇x 1 , see e.g [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF][START_REF] Condat | Discrete total variation: New definition and minimization[END_REF]; it yields images with sharp edges, but the textures are still over-smoothed, there are staircasing effects and the pixel values tend to be clustered in piecewise constant areas. To overcome these drawbacks, one can penalize higher order derivatives [START_REF] Lefkimmiatis | Hessian Schatten-norm regularization for linear inverse problems[END_REF][START_REF] Bredies | Total generalized variation[END_REF] or make use of nonlocal penalties [START_REF] Peyré | Non-local regularization of inverse problems[END_REF][START_REF] Condat | Semi-local total variation for regularization of inverse problems[END_REF][START_REF] Chierchia | A nonlocal structure tensor-based approach for multicomponent image recovery problems[END_REF]. Another approach, which is at the heart of the recent paradigm of sparse recovery [START_REF] Starck | Sparse image and signal processing: wavelets, curvelets, morphological diversity[END_REF][START_REF] Elad | Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing[END_REF] and compressed sensing [START_REF] Donoho | Compressed sensing[END_REF][START_REF] Strohmer | Measure what should be measured: progress and challenges in compressive sensing[END_REF], is to choose R to favor some notion of low complexity. Indeed, many phenomena, when observed by instruments, yield data living in high dimensional spaces, but inherently governed by a small number of degrees of freedom. One early choice was to set R as the 1 norm of wavelet coefficients of the image. But the signals encountered in applications like radar, array processing, communication, seismology, or remote sensing, are usually specified by parameters in a continuous domain, from which they depend nonlinearly. So, modern sampling theory has widened its scope to a broader class of signals, with so-called finite rate of innovation, i.e. ruled by parcimonious models [START_REF] Maravic | Sampling and reconstruction of signals with finite rate of innovation in the presence of noise[END_REF][START_REF] Dragotti | Sampling moments and reconstructing signals of finite rate of innovation: Shannon meets Strang-Fix[END_REF][START_REF] Blu | Sparse sampling of signal innovations[END_REF][START_REF] Urigiien | Sampling at the rate of innovation: Theory and applications[END_REF]. This encompasses reconstruction of pulses from lowpass measurements [START_REF] Condat | Cadzow denoising upgraded: A new projection method for the recovery of Dirac pulses from noisy linear measurements[END_REF] and spectral estimation, which is the reconstruction of sinusoids from point sam-ples [START_REF] Stoica | List of references on spectral line analysis[END_REF][START_REF] Stoica | Spectral Analysis of Signals[END_REF], with many applications [START_REF] Chan | A parameter estimation approach to estimation of frequencies of sinusoids[END_REF][START_REF] Johnson | The application of spectral estimation methods to bearing estimation problems[END_REF][START_REF] Tufts | Estimation of frequencies of multiple sinusoids: Making linear prediction perform like maximum likelihood[END_REF][START_REF] Bresler | Exact maximum likelihood parameter estimation of superimposed exponential signals in noise[END_REF][START_REF] Stoica | Maximum likelihood estimation of the parameters of multiple sinusoids from noisy measurements[END_REF][START_REF] Macleod | Fast nearly ML estimation of the parameters of real or complex single tones or resolved multiple tones[END_REF][START_REF] Quinn | The estimation and tracking of frequency[END_REF][START_REF] So | Linear prediction approach for efficient frequency estimation of multiple real sinusoids: algorithms and analyses[END_REF][START_REF] Condat | A 2-D spectral analysis method to estimate the modulation parameters in structured illumination microscopy[END_REF]. The knowledge of the kind of elements we want to promote in the image makes it possible to estimate them from coarse-scale measurements, even with infinite precision if there is no noise. Methods achieving this goal are qualified as super-resolution methods, because they uncover fine scale information, which was lost in the data, beyond the Rayleigh or Nyquist resolution limit of the acquisition system [START_REF] Fernandez-Granda | Super-resolution and compressed sensing[END_REF][START_REF] Candès | Towards a mathematical theory of super-resolution[END_REF]. However, in this context, maximum likelihood estimation amounts to structured low rank approximation, which forms nonconvex and very difficult, even NP-hard in general, problems [START_REF] Markovsky | Low rank approximation: algorithms, implementation, applications[END_REF][START_REF] Condat | Cadzow denoising upgraded: A new projection method for the recovery of Dirac pulses from noisy linear measurements[END_REF]. An elegant and unifying formulation, which yields convex problems, is based on the atomic norm [START_REF] Bhaskar | Atomic norm denoising with applications to line spectral estimation[END_REF][START_REF] Condat | Atomic norm minimization for decomposition into complex exponentials[END_REF]. We place ourselves in this general framework of atomic norm minimization: the sought-after image x is supposed to be a sparse positive combination of the elements of an infinite dictionary A, indexed by continuously varying parameters. Then, one can choose R as the atomic norm x A of the image x, which can be viewed as the 1 norm of the coefficients, when the image is expressed in terms of the unit-norm elements of A, called atoms:

(1)

x A = inf {t > 0 : x ∈ t conv(A)} , where conv(A) is the convex hull of the atoms. In this paper, we consider the setting, which is new to our knowledge, where the atoms are lines. Expressed in the Fourier domain, these atoms can be characterized with respect to their rows and columns and the problem can be reduced to a dictionary of 1-D complex exponential samples, indexed by their frequency and phase, whose atomic norms can be computed via semidefinite programming [START_REF] Vandenberghe | Semidefinite programming[END_REF]. This formulation makes it possible to derive a convex optimization problem under constraints, solved by mean of a primal-dual splitting algorithm [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]. Then, applying a Prony-like method [START_REF] Rahman | Total least squares approach for frequency estimation using linear prediction[END_REF] to the solution of the algorithm allows us to extract the parameters of the lines. This approach estimates the lines with high accuracy, whereas the Hough [START_REF] Hough | Method and means for recognizing complex patterns[END_REF][START_REF] Illingworth | A survey of the Hough transform[END_REF][START_REF] Mukhopadhyay | A survey of Hough transform[END_REF] and the Radon [START_REF] Radon | 1.1 über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeitenestimmung[END_REF][START_REF] Deans | Hough transform from the Radon transform[END_REF][START_REF] Murphy | Linear feature detection and enhancement in noisy images via the Radon transform[END_REF] transforms fail, due to their discrete nature. Our motivation stems from the frequent presence in biomedical images, e.g. in microscopy, of elongated structures like filaments, neurons, veins, which are deteriorated when reconstructed with classical penalties. Some related works are dedicated to the recovery of curve-like singularities, by variational methods [START_REF] Aubert | Analysis of a new variational model to restore point-like and curve-like singularities in imaging[END_REF][START_REF] Yu | Geometric mode decomposition[END_REF], Riesz-based models [START_REF] Marchant | Using super-resolution methods to solve a novel multisinusoidal signal model[END_REF][START_REF] Marchant | A sinusoidal image model derived from the circular harmonic vector[END_REF] or so-called finite rate of innovation methods [START_REF] Maravic | Sampling and reconstruction of signals with finite rate of innovation in the presence of noise[END_REF][START_REF] Shukla | Sampling schemes for multidimensional signals with finite rate of innovation[END_REF][START_REF] Chen | 2D finite rate of innovation reconstruction method for step edge and polygon signals in the presence of noise[END_REF][START_REF] Pan | Sampling curves with finite rate of innovation[END_REF][START_REF] Fatemi | Sampling and reconstruction of shapes with algebraic boundaries[END_REF][START_REF] Ongie | Off-the-grid recovery of piecewise constant images from few Fourier samples[END_REF][START_REF] Pan | Sampling curves with finite rate of innovation[END_REF]. The originality of our method is that we reduce the minimization over an infinite dictionary of lines to a semidefinite programming problem, taking advantage of the line structure in both directions of the grid. Although there are works in the same vein to recover 2-D point sources, or equivalently to estimate the parameters of 2-D exponentials [START_REF] Potts | Parameter estimation for exponential sums by approximate Prony method[END_REF][START_REF] Plonka | How many Fourier samples are needed for real function reconstruction?[END_REF][START_REF] Cuyt | Multivariate exponential analysis from the minimal number of samples[END_REF][START_REF] Dossal | Sampling the Fourier transform along radial lines[END_REF], applying similar principles to the estimation of lines is not straightforward and is new, to our knowledge.

The paper is organized as follows. The model is described in section 2, the framework of atomic norm minimization underlying the super-resolution principle is introduced in section 3, the algorithms we derive are presented in section 4. Then a Prony-like method is developed in section 5, as a way to perform spectral estimation of the line parameters. Section 6 gives an overview of related works. Finally some experimental results are shown in section 7. Part of this work has been published in a conference paper [START_REF] Polisano | Convex super-resolution detection of lines in images[END_REF]. In the present paper, we added mathematical developments, another algorithm, a new estimation procedure for the line parameters, an extension to the whole range of line angles without restriction any more, an application to inpainting problems and several numerical experiments.

Notations. Vectors, e.g. z = (z 0 , . . . , z N -1 ) and matrices have component indexes starting at zero. The Hilbert space of complex matrices of size M × N is denoted by M M,N . The entry in the (k 1 + 1)-th row and (k 2 + 1)-th column of a matrix M is referred to M[k 1 , k 2 ] or M k1,k2 . Due to required multiple subscripts, we often adopt the matlab™ notation M[k 1 , :] (resp. M[:, k 2 ]) for referring to the (k 1 +1)-th row (resp. (k 2 +1)-th column) of the matrix M. To extract a submatrix, we also use the notation M[p 1 : q 1 , p 2 : q 2 ]. The nuclear norm of the matrix M, denoted by M * = i σ i (M), is the sum of its singular values σ i (M). The operation M ⊗ N denotes the Kronecker product between a m × n matrix M and a p × q matrix N, whose result is the mp × nq matrix (2)

M ⊗ N =    M 0,0 N • • • M 0,n-1 N . . . . . . . . . M m-1,0 N • • • M m-1,n-1 N    .
The ground truths for distributions or matrices to recover are denoted by a sharp • , while the estimated variables are denoted by a tilde •. Continuous (resp. discrete) 1-D Fourier transform of functions or distributions (resp. vectors) is denoted by a hat •. For 2-D functions, distributions, or matrices, the notation f always refers to the horizontal Fourier transform of f , while the 2-D Fourier transform is denoted by Ff . When it is necessary to distinguish them, the horizontal and vertical Fourier transform are denoted by F 1 f and F 2 f , respectively. The Radon transform of a function f is denoted by Rf . The projection operator onto a set C is denoted by P C and the proximity operator of a closed proper convex function f : R N → R ∪ {+∞} (see [START_REF] Boyd | Convex optimization[END_REF][START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Parikh | Proximal algorithms[END_REF]) is defined by

(3) prox f (x) = arg min y f (y) + 1 2 x -y 2 2 .
2. An image model of blurred lines. Our aim is to restore lines from an observed image y = b + , which is made of a blurred image b containing lines and corrupted by some noise and then to estimate the parameters-angle, offset, amplitude-of the lines. In this section, we formulate what we precisely mean by an image containing lines. In short, b is a sum of perfect lines which have been blurred and then sampled. Both processes are detailed in the following.

2.1. The ideal continuous model and the objectives. We place ourselves in the quotient space P = R/(W Z) × R, corresponding to the 2-D plane with horizontal W -periodicity, for some integer W 1. To simplify the notations, we suppose that W is odd and we set M = (W -1)/2.

A perfect line, with angle θ ∈ (-π/2, π/2] with respect to the vertical axis, amplitude α > 0 and offset η ∈ R from the origin on the horizontal axis, can be defined as a tempered distribution, which maps a function ψ in the Schwartz class S(R 2 ) to its integral along the geometric line L = {(t 1 , t 2 ) ∈ P : (t 1 -η) cos θ + t 2 sin θ = 0} ; that is, to α L ψ(t 1 , t 2 ) dt 1 dt 2 . Thus, by a slight abuse of notation, we can write the perfect line as:

(t 1 , t 2 ) ∈ P → αδ (t 1 -η) cos θ + t 2 sin θ ,
where δ is the (1-D) Dirac distribution. For more details about multi-dimensional distributions over curves or surfaces, we refer the interested reader to [START_REF] Gelfand | Properties and operations[END_REF][START_REF]Impulse functions over curves and surfaces and their applications to diffraction[END_REF]. We define the distribution x , which is a sum of K different such perfect lines, for some integer K 1, as [START_REF] Bhaskar | Atomic norm denoising with applications to line spectral estimation[END_REF] x

t 2 t 1 θ η (a) t 2 t 1 n 2 n 1 (b) (c)
: (t 1 , t 2 ) ∈ P → K k=1 α k δ (t 1 -η k ) cos θ k + t 2 sin θ k .
Figure 2a illustrates the line parameters and Figure 2b the convention we use for representing images.

Hypothesis 1. At this time, we suppose that the lines are rather vertical; that is, for every

k = 1, . . . , K, θ k ∈ (-π/4, π/4].
This hypothesis is made, because the rows and columns of the image will be processed differently. We will proceeed in the Fourier domain, by applying the discrete 1-D Fourier transform on every row of the image; we insist on the fact that we do not consider the usual 2-D Fourier transform. This setting may appear restrictive, but we show in subsection 4.3 that, in fact, we are able to deal with the general case.

2.2.

A blur model for an exact sampling process. The observed image b of size W × H in Figure 2c is obtained by convolution of the distribution x with a blur function φ, followed by sampling with unit step:

(5) b [n 1 , n 2 ] = (x * φ)(n 1 , n 2 ), ∀n 1 = 0, . . . , W -1, n 2 = 0, . . . , H -1 .
The point spread function (PSF) φ is assumed to be separable; that is, it can be written φ(t 1 , t 2 ) = ϕ 1 (t 1 )ϕ 2 (t 2 ). So, the function b = x * φ can be obtained by a first horizontal convolution with ϕ 1 and then a second vertical convolution with ϕ 2 . We can show (see Appendix A) that we get the function

(6) b = x * φ : (t 1 , t 2 ) ∈ P → K k=1 α k ψ k (t 1 -η k ) cos θ k + t 2 sin θ k , where (7)
ψ k = 1 cos θ k ϕ 1 • cos θ k * 1 sin θ k ϕ 2 • sin θ k , if θ k = 0 and ψ k = ϕ 1 else.
Remark 1. We can notice that (6) can also be interpreted as follows: every line has undergone a 1-D convolution with ψ k in the direction transverse to it. We can also notice that if ϕ 1 and ϕ 2 are Gaussian functions and have same variance κ 2 , it follows from (7) that ψ k has variance κ 2 cos 2 θ + sin 2 θ = κ 2 as well.

Assumptions. We assume that ϕ 1 and ϕ 2 have the following properties:

(i) ϕ 1 ∈ L 1 (0, W ) is W -periodic, bounded, such that 1 W W 0 ϕ 1 (t 1 ) dt 1 = 1 and bandlimited; that is, its Fourier coefficients c m (ϕ 1 ) = 1 W W 0 ϕ 1 (t 1 )e -j2πmt1/W dt 1
are zero for every m ∈ Z with |m| (W + 1)/2 = M + 1. Then, the discrete filter

(8) g[n] = ϕ 1 (n) n∈Z has discrete Fourier coefficients g[m] = 1 W W -1 n=0 g[n]e -j 2πmn W = c m (ϕ 1 ) . (ii) ϕ 2 ∈ L 1 (R) is such that R ϕ 2 (t 2
) dt 2 = 1 and bandlimited so that the discrete filter [START_REF] Bredies | Total generalized variation[END_REF] h

[n] = ϕ 2 (n) n∈Z
has compact support of length 2S + 1, for some S ∈ N; that is,

h[n] = 0, if |n| S + 1 .
An example of three blurred lines is depicted in Figure 1a. We insist on the fact that no discrete approximation is made during the sampling process leading from the continuous to the discrete formulation, due to the Assumptions (i), (ii) and Hypothesis 1. These important guarantees, for the purpose of super-resolution, are detailed in Appendix A.

Remark 2. These assumptions are convenient and reasonable. Indeed, in practice, we can always approximate a PSF by a bandlimited function. That said, we also propose in Appendix A a weaker version of Assumption (ii), denoted by (ii'), which relax the bandlimited restriction. Regarding the periodicity in Assumption (i), this hypothesis can be circumvented by considering an acquisition process yielding a twice larger image, in order to "periodize" the observed image. To achieve this, let us also denote by x the discrete image obtained by sampling x = F 1 x , which is the horizontal Fourier transform, in the sense of distributions, of the ideal model x (4) made up of 1-D Dirac distributions. Therefore, x is composed of a sampled sum of exponentials:

x [m, n 2 ] = K k=1 α k cos θ k e j2π(n2 tan θ k -η k ) m W , (10)
m = -M, . . . , M, n 2 = -S, . . . , H -1 + S . Remark 3. x is such that x [-m, n 2 ] = x [m, n 2 ] *
, where • * denotes the complex conjugation, so we can only deal with the right part x [0 : M, : ] and we can note that the row corresponding to m = 0 is real and equal to K k=1 α k cos θ k . We consider in the following the Fourier image x [m, n 2 ] of size (M + 1) × H S , with H S = H + 2S, due to the addition of S pixels beyond the boundaries, allowing the convolution with filter h determined hereafter.

Our goal will be to reconstruct x from b , that is from its observations through a known degradation operator A, which we characterize in the following proposition: Proposition 1. Let A be the operator which for all m = 0, . . . , M multiplies each row vector x [m, :] of x by the corresponding Fourier coefficient g[m] and convolves it with the filter h = (h -S , . . . , h 0 , . . . , h S ). Then, we have the relation

(11) A x = b ,
which can be alternatively rewrite on each row m = 0, . . . , M as:

(12) g[m]( x [m, :] * h) = b [m, :] . Proof. See in Appendix B.
Remark 4. From a matrix point of view, the operator A corresponds to a left and right matrix multiplication with the matrices G of size (M + 1) × (M + 1) and the transpose of Ȟ of size H × H S defined by:

(13) G = diag( g 0 , . . . , g M ), Ȟ =        h -S • • • h S 0 0 • • • 0 0 h -S • • • h S 0 • • • 0 . . . • • • . . . . . . . . . • • • . . . 0 • • • 0 h -S • • • h S 0 0 • • • 0 0 h -S • • • h S        ; that is, (14)
A x = G x ȞT .
Finally, the image b of the blurred lines is corrupted by noise, so that we observe the degraded image [START_REF] Carathéodory | Über den Zusammenhang der Extremen von harmonischen Funktionen mit ihren Koeffizienten und über den Picard-Landau'schen Satz[END_REF]

y = b + , ∼ N (0, ζ) ,
with ζ the noise level. Thus, the problem consists from ( 15) and [START_REF] Briani | VEXPA: Validated EXPonential Analysis through regular sub-sampling[END_REF] in recovering x from y, which takes the form of a ill-posed linear inverse problem. Then, we will need to express a convex optimization problem, under constraints exploiting the sparse structure of the signal we are looking for, namely that it is a combination of lines. The super-resolution process consists of recovering the high frequency content (lost because of the blur operator) from the degraded image y, which can be viewed as a spectral extrapolation process. Then, we aim at recovering the parameters (θ k , η k , α k ) of these lines, as a post-process, after the reconstruction. This procedure will be decomposed as:

1. First, solve a convex optimization problem of the form:

(16) Minimize y -A x , under the constraint that x is made of lines ;

that is, to go to the bottom line of the diagram in Figure 15 from y to x . 2. Second, perform a Prony-like method on x in order to estimate the K parameters (θ k , η k , α k ).

These two steps are summarized in Figure 3. Note that this work also covers the case where a mask is applied; that is, it can encompass inpainting problems. In the next section, we present the framework of atomic norm, from which the optimization problem will be derived.

Super-resolution detection of lines.

3.1. Atomic norm and semidefinite characterizations. Consider a complex signal z ∈ C N represented as a K-sparse mixture of atoms from the set

A = a(ω) ∈ C N : ω ∈ Ω , that is, z = K k=1 c k a(ω k ), c k 0, ω k ∈ Ω .
We consider atoms a(ω) ∈ C N that are continuously indexed in the dictionary A by the parameter ω in a compact set Ω. The atomic norm, first introduced in [START_REF] Chandrasekaran | The convex geometry of linear inverse problems[END_REF], is defined as

z A = inf {t > 0 : z ∈ t conv(A)} ,
where conv(A) denotes the convex hull of a general atomic set A, enforcing sparsity. Chandrasekaran et al. [START_REF] Chandrasekaran | The convex geometry of linear inverse problems[END_REF] argue that the atomic norm is the best convex heuristic for underdetermined, structured linear inverse problems, which generalizes the 1 norm for sparse recovery and the nuclear norm for low-rank matrix completion. In our problem (4) the atoms are somehow lines, so one can considered in the Fourier domain the dictionary A 2D indexed by the angle and the offset; that is, composed of the 2-D exponential atoms of size W × H S :

a 2D (θ k , η k ) = 1 cos θ k e j2π(n2 tan θ k -η k )m/W , m = -M, . . . , M, n 2 = -S, . . . , H -1 + S ,
as illustrated in Figure 4. The problem is that there is no closed-form expression for the atomic norm in these 2-D dictionaries, to our knowledge. However, in the case of 1-D complex exponentials, there is a way to compute the atomic norm via semidefinite programming. So, we reformulate the problem using the simplified 1-D case. From now on, we consider the dictionaries

A = a(f, φ) ∈ C |I| , f ∈ [0, 1], φ ∈ [0, 2π) , (17)
A 0 = a(f ) ∈ C |I| , f ∈ [0, 1] , (18)
in which the atoms are the vectors of components [a(f, φ)] i = e j(2πf i+φ) , i ∈ I and

= α 1 + α 2 + α 3 = α 1 + α 2 + α 3
Fig. 4. Illustration with a signal made of a weighted combination of three lines (in gray). In the Fourier domain, we have through the sampling process the same kind of combination, but with 2-D exponential atoms a 2D (θ k , η k ). In both cases, the weights α i are the amplitudes of the lines.

[a(f )] i = [a(f, 0)] i = e j2πf i , i ∈ I. The atomic norm writes: [START_REF] Chandrasekaran | The convex geometry of linear inverse problems[END_REF] z A = inf

c k 0 f k ∈[0,1) φ k ∈[0,2π) k c k : z = k c k a(f k , φ k ) .
Theorem 2 (Caratheodory). A vector z = (z * N -1 , . . . , z * 1 , z 0 , z 1 , . . . , z N -1 ) of length 2N -1, with z 0 ∈ R, is a positive combination of K N atoms a(f k ) if and only if T N (z + ) 0 and is of rank K, where z + = (z 0 , . . . , z N -1 ) is of length N , T N : C N → T N is the Toeplitz operator (20)

T N : z + = (z 0 , . . . , z N -1 ) →      z 0 z * 1 • • • z * N -1 z 1 z 0 • • • z * N -2 . . . . . . . . . . . . z N -1 z N -2 • • • z 0     
, with T N the Hilbert subspace of M N composed of Hermitian Toeplitz matrices and 0 denotes positive semidefiniteness. Moreover, this decomposition is unique, if K < N .

Proof. See references [START_REF] Carathéodory | Über den Variabilitätsbereich der Fourier'schen Konstanten von positiven harmonischen Funktionen[END_REF][START_REF] Carathéodory | Über den Zusammenhang der Extremen von harmonischen Funktionen mit ihren Koeffizienten und über den Picard-Landau'schen Satz[END_REF][START_REF] Toeplitz | Zur Theorie der quadratischen und bilinearen Formen von unendlichvielen Veränderlichen[END_REF][START_REF] Polisano | Modélisation de textures anisotropes par la transformée en ondelettes monogènes et super-résolution de lignes 2-D[END_REF].

We are now in position to characterize the atomic norm:

Proposition 3. The atomic norm z A can be characterized by the following semidefinite program:

(21) z A = min q∈C N ,q0 0 q 0 : T N (z, q) = T N (q) z z * q 0 0 ,
where q 0 is the first component of vector q = (q 0 , . . . , q N -1

) ∈ R + × C N -1 and z * = z T .
Proof. This result is an improvement of [82, Proposition II.1] and the proof is given in Appendix C.

Since the matrix T N (z, q) in ( 21) is Hermitian and positive semidefinite, its eigenvalues (λ i ) 0 i N are positive reals. So, q 0 = 1 N +1 tr(T N (z, q)) = 1 N +1 N i=0 λ i is real and positive.

3.2. Properties of the model x with respect to the atomic norm. In Fourier domain, the discrete image x given by [START_REF] Bresler | Exact maximum likelihood parameter estimation of superimposed exponential signals in noise[END_REF] can be viewed as a sum of atoms: regarding the columns l n2 of the matrix x , with I = {0, . . . , H S -1}, we have:

(22) l n2 = x [:, n 2 ] = K k=1 c k a(f n2,k ) ,
and regarding the rows t m , with I = {-M, . . . , M }, we have:

(23) t m = x [m, :] = K k=1 c k a(f m,k , φ m,k ) T ,
where

c k = α k cos θ k , f n2,k = n 2 tan θ k -η k W , (24)
φ m,k = - 2πmη k W , f m,k = m tan θ k W , d m,k = c k e jφ m,k , e m,k = e jφ m,k .
The vectors l n2 of size W = 2M + 1 are positive combinations of K atoms a(f n2,k ), with K M , since we can reasonably assume that the number of lines K is smaller than half the number of pixels M . Thus, Theorem 2 ensures that the decomposition ( 22) is unique, hence, following ( 19):

(25)

l n2 A = K k=1 c k = x [0, n 2 ], ∀n 2 = 0, . . . , H S -1 .
By contrast, since the coefficients d m,k are complex, Theorem 2 no longer holds and we simply have from Proposition 3:

(26) t m A K k=1 c k , ∀m = -M, . . . , M .
Let us take a closer look at the case of one line; that is, K = 1, characterized by parameters (θ, η, α). We recall by ( 10) that x can be written as:

x [m, n 2 ] = c 1 e j2π((f1-f0)n2+f0)m , c 1 = α cos θ , f 0 = - η W , f 1 = tan θ -η W .
Let z = (z 0 , . . . , z N -1 ) be a complex vector, whose elements z i are rearranged in a Toeplitz matrix P K (z) of size (N -K) × (K + 1) and rank K as follows:

P K (z) =    z K • • • z 0 . . . . . . . . . z N -1 • • • z N -K-1    .
We get the following characterization of one line in Fourier domain:

Proposition 4. An image x is of the form x[m, n] = c 1 e j2π((f1-f0)n+f0)m)
if and only if the columns l n and rows t m of x are such that T M (l n ) is positive semidefinite and of rank one, P 1 (t m ) is of rank one and x[0, n] = x[0, 0] for all m and n.

Proof. See Appendix D. Besides, with D = diag(c 1 , . . . , c K ) and V n2 = a(f n2,1 ) • • • a(f n2,K ) , we can remark that T M (l n2 ) = K k=1 c k T M (a(f n2,k )) = K k=1 c k a(f n2,k )a(f n2,k ) * = V n2 DV * n2 ,
where • * denotes the hermitian conjugate. Since the Toeplitz matrices T M (a(f n2,k )) only contain ones on their main diagonal, then

1 M tr(T M (l n2 )) = K k=1 c k .
Moreover, the trace of a positive semidefinite matrix is equal to its nuclear norm, hence:

T M (l n2 ) * = tr(T M (l n2 )) = M K k=1 c k = M l n2 A .
The nuclear norm of a matrix, which corresponds to the sum of its singular values, is often used as a convex approximation of the rank of this matrix [START_REF] Recht | Necessary and sufficient conditions for success of the nuclear norm heuristic for rank minimization[END_REF][START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF]. Consequently, in the following, we consider a convex relaxation of the line characterization given in Proposition 4, in which the rank constraint on T M (l n2 ) is replaced by an atomic norm constraint on l n2 A . Since the minimum value to achieve is c = K k=1 c k and since the atomic norm lies on the first row x[0, n 2 ] = x[0, 0], we impose the constraint x[0, n 2 ] = x[0, 0] c . We do the same for the rows.

Remark 5. We showed that if every row and column of the horizontal Fourier transform of the image satisfy assumptions of the Proposition 4, then we have the horizontal Fourier transform of an image containing one line. So, Proposition 4 supports our strategy of dealing with the rows and columns of the image, to reformulate the 2-D problem as a combination of 1-D problems. Proposition 4 shows that, in the case of one line, we do not lose anything by this process; in other words, there is no other image than the image containing one line, whose rows and columns have the prescribed form. This is an indication (not a proof ) that the upcoming SDP formulation (27)-( 28) is tight: we will not promote structures other than lines in the image. We could not derive a similar characterization for K 2 lines, which is a difficult task, but the philosophy remains: the aim is to minimize atomic norms of rows and columns simultaneously, so that the solution will be composed of sparse sums of exponentials in both directions. The convex optimization problem exploiting this strategy is presented in the next section.

4. Minimization problem with atomic norm regularization. Given the operator A : X → Y defined in ( 14) using the filters ( 8)-( 9) and y the Fourier version of the degraded image observed (15), we are looking for an image x ∈ X which minimizes A xy Y , for the norm derived from the inner product [START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF] and whose rows and columns satisfy properties [START_REF] Condat | Atomic norm minimization for decomposition into complex exponentials[END_REF] and [START_REF] Condat | Discrete total variation: New definition and minimization[END_REF]. We fix a constant c c . Consequently, the following optimization problem provides an estimator of x defined in (10): [START_REF] Condat | A 2-D spectral analysis method to estimate the modulation parameters in structured illumination microscopy[END_REF] x ∈ arg min

( x,q)∈X ×Q 1 2 A x -y 2 Y , s.t              x[0, n 2 ] = x[0, 0] c , (28a) q[m, 0] c , (28b) T M +1 ( x[:, n 2 ]) 0 , (28c) T H S ( x[m, :], q[m, :]) 0 , ( 28d
)
∀n 2 = 0, ..., H S -1, ∀m = 1, ..., M ,
where the real Hilbert spaces (X ,

•, • X ), (Y, •, • Y ) and (Q, •, • Q ) are respectively defined by (29) X = x ∈ C (M +1)×H S : Im( x[0, :]) = 0 ,
endowed with the following inner product:

(30) x 1 , x 2 X = H S -1 n2=0 x 1 [0, n 2 ] x 2 [0, n 2 ] + 2Re M m=1 H S -1 n2=0 x 1 [m, n 2 ] x 2 [m, n 2 ] * ; (31) Y = y ∈ C (M +1)×H : Im( y[0, :]) = 0 ,
which is equivalent to [START_REF] Cuyt | Multivariate exponential analysis from the minimal number of samples[END_REF] for S = 0, endowed with the inner product:

(32) y 1 , y 2 Y = H-1 n2=0 y 1 [0, n 2 ] y 2 [0, n 2 ] + 2Re M m=1 H-1 n2=0 y 1 [m, n 2 ] y 2 [m, n 2 ] * ; (33) Q = q ∈ C (M +1
)×H S : Im(q[:, 0]) = 0 , endowed with the following inner product:

(34) q 1 , q 2 Q = M m=0 q 1 [m, 0]q 2 [m, 0] + 2Re H S -1 n2=1 M m=0 q 1 [m, n 2 ]q 2 [m, n 2 ] * .
Remark 6. Since we deal with rows and columns of x ∈ X and rows of q ∈ Q, we also mention the real Hilbert spaces they belong to, which are respectively denoted by

X l ⊂ R × C M , X t ⊂ C H S and Q t ⊂ R × C H S -1 endowed with the following inner products: l 1 , l 2 X l = l 1 [0]l 2 [0] + 2Re M m=1 l 1 [m]l 2 [m] * , (35) t 1 , t 2 Xt = 2Re H S n2=1 t 1 [n 2 ]t 2 [n 2 ] * , ( 36
)
q 1 , q 2 Qt = q 1 [0]q 2 [0] + 2Re H S -1 n2=1 q 1 [n 2 ]q 2 [n 2 ] * . ( 37
)
The operators T M +1 : X l → T M +1 and T H S : X t × Q t → T H S +1 are defined respectively on columns and rows of x ∈ X and q ∈ Q, endowed respectively with the inner products ( 35) and ( 36)- [START_REF] Dragotti | Sampling moments and reconstructing signals of finite rate of innovation: Shannon meets Strang-Fix[END_REF], to Hermitian Toeplitz matrices of dimension M + 1 and H S + 1 respectively, whose spaces are denoted by T M +1 and T H S +1 endowed with the classical inner product on complex matrices:

(38) M, N M = i,j M * ij N ij ,
and corresponding Frobenius norm

(39) M F =   i,j |M ij | 2   1/2
.

The expressions of the operators T N and T N are given respectively in ( 20) and ( 21).

Remark 7. This optimization problem could be rewritten in a regularized form involving a parameter λ to tune, what is not any better than the tuning parameter c, which has the advantage of having a physical meaning, related to the line intensities. Moreover, we have a simple estimator of c noticing that

b [0, n 2 ] = g[0] x [0, :] * h [n 2 ] = i h[i] g[0]c = c , ∀n 2 = 1, . . . , H . Then, since y[0, n 2 ] = b [0, n 2 ]
+ with E( ) = 0, one can get an estimation of the parameter c by averaging the first row:

c ≡ 1 H H-1 n2=0 y[0, n 2 ] ≈ c .
We keep this constrained formulation and write it in a more suitable way as follows.

Let H = X × Q be the Hilbert space in which the variable X = ( x, q) lies, endowed with the following inner product:

(40) ( x 1 , q 1 ), ( x 2 , q 2 ) H = x 1 , x 2 X + q 1 , q 2 Q . Let us define L (1)
m : H → T H S +1 and L

(2)

n2 : H → T M +1 by L (1) m (X) = T H S ( x[m, :], q[m, :]) , (41)
L (2) n2 (X) = T M +1 ( x[:, n 2 ]) . ( 42
)
We denote by ι C the indicator function of a convex set C, defined by

ι C : x → 0 ifx ∈ C +∞ ifx / ∈ C .
We denote by C the cone of positive semidefinite matrices and we introduce the set B ⊂ H corresponding to the boundary constraints:

(43) B = ( x, q) ∈ H : x[0, n 2 ] = x[0, 0] c, q[m, 0] c .
Then the optimization problem [START_REF] Condat | A 2-D spectral analysis method to estimate the modulation parameters in structured illumination microscopy[END_REF] under constraints (28) can be rewritten as follows: (44)

X = arg min X=( x,q)∈H 1 2 A x -y 2 Y + ι B (X) + M m=1 ι C (L (1) m (X)) + H S -1 n2=0 ι C (L (2) n2 (X)) .
We now propose two different algorithms to solve this convex optimization problem. The first one in subsection 4.1 is more general and can be applied to the extended setting (with no restriction on the line angles) and inpainting problems, presented in subsection 4.3 and subsection 4.4 respectively, whereas the second cannot, since it is difficult to compute the proximity operator associated to the data fidelity term. However, the second algorithm detailed in subsection 4.2, happens to be faster and so will be used for most of the numerical experiments exposed in section 7.

4.1. First algorithm design. The optimization problem [START_REF] Gelfand | Properties and operations[END_REF] can be cast as a minimization problem, involving smooth, proximable and linear composite terms [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]:

(45) X = arg min X∈H F (X) + G(X) + Q-1 i=0 H i (L i (X)) , with F (X) = 1 2 A x -y 2 Y , X = ( x, q), G = ι B , which is proximable, Q = M + H S linear composite terms where H i = ι C , L i = L (2) i when 0 i H S -1 and L i = L (1) i-H S +1 when H S i H S + M -1. We define Hx = Q-1 i=0 H i x i , L (1) (X) = (L (1) 1 (X), . . . , L (1)
M (X)) and L (2) (X) = (L (2) 0 (X), . . . , L (2)
H S -1 (X)). L = (L (1) , L (2) ) is the linear operator such that the composite terms can be written H • L, where • denotes composition. We define an inner product on the range of L (1) , L (2) and L, which are product spaces, as the sum of the inner products defined on these spaces (similarly to [START_REF] Eftekhari | A bridge between past and present: Exchange and conditional gradient methods are equivalent[END_REF]). We denote the corresponding norms by • (1) , • (2) and • (1,2) . We define the following operator norms:

A = sup x∈X A x Y x X , (46)
L i = sup X∈H L i (X) F X H , (47)
L (j) = sup X∈H L (j) (X) (j) X H , j ∈ {1, 2} , (48)
L = sup X∈H L(X) (1,2) X H . ( 49
)
We now establish some properties of these functions, operators and norms.

Lemma 5. The norm (46) of the operator A defined in [START_REF] Carathéodory | Über den Variabilitätsbereich der Fourier'schen Konstanten von positiven harmonischen Funktionen[END_REF] is given by

(50) A = g ∞ h ∞ .
Proof. If we denote by x k the (k+1)-th row of x, then by definition the operator A maps the (k+1)-th row of x to g[k]( x k * h). Let us calculate the norm of this operator. By considering in Fourier the norm operator f → f * h we have the inequality

x k * h 2 h ∞ x k 2 .
Thus, with the norm derived from (30) we get

A x 2 Y = | g 0 | 2 x 0 * h 2 2 + 2| g 1 | 2 x 1 * h 2 2 + • • • + 2| g M | 2 x M * h 2 2 , (51)
g 2 ∞ h 2 ∞ ( x 0 2 2 + 2 x 1 2 2 + • • • + 2 x M 2 2 ) , g 2 ∞ h 2 ∞ x 2 X .
Since the filter h is lowpass, the equality is attained for an image x whose all rows are zero, except one which is constant (and nonzero), of index m 0 , where g m0 corresponds to the maximum g ∞ , which proves the result [START_REF] Lefkimmiatis | Hessian Schatten-norm regularization for linear inverse problems[END_REF].

Lemma 6. The adjoint operator of the operator A is denoted by A * and defined such that A x , z Y = x , A * z X . Its matrix expression is given by

(52) A * z = G * z Ȟ .
Proof. From the definition ( 14), we have the matrix product A x = G x ȞT . Then we use the fact that for any matrix M, we have M x 1 , x 2 M = x 1 , M * x 2 M and we remark that

x 1 , x 2 X = x 1 , x 2 M + x 1 , x 2 * M -x 1 [0, :] , x 2 [0, :] R H S .
This yields the claimed result.

Then, we have the following proposition:

Proposition 7. For X = ( x, q), the gradient of F (X) = 1 2 A x -y 2 Y is ∇F (X) = (A * (A x -y), 0) T , which is Lipschitz-continuous with Lipschitz constant β = g 2 ∞ h 2 ∞ . Proof. See in Appendix G.
We now give in the next proposition the adjoint operators of:

T M +1 : (X l , • , • X l ) → (T M +1 , • , • M ) , (53)
T H S : (X t × Q t , • , • Xt×Qt ) → (T H S +1 , • , • M ) , (54)
where the inner products are defined in [START_REF] Donoho | Compressed sensing[END_REF], [START_REF] Dossal | Sampling the Fourier transform along radial lines[END_REF], [START_REF] Dragotti | Sampling moments and reconstructing signals of finite rate of innovation: Shannon meets Strang-Fix[END_REF] and [START_REF] Eftekhari | Stable super-resolution of images[END_REF]. Proposition 8. For M (1) ∈ T H S +1 and M (2) ∈ T M +1 the adjoint operators of (53) and (54) applied to M (1) and M (2) give the vectors

z 2 = T * M +1 M (2) ∈ R × C M , (z 1 , q 1 ) = T * H S M (1) ∈ C H S × (R × C H S -1 ) ,
respectively, whose components are:

z 2 [k] = M -k l=0 M (2)
l+k,l , ∀k = 0, . . . , M ,

z 1 [k] = M (1)
H S +1,k , q 1 [k] = H S -1-k l=0 M (1) l+k,l + δ k M (1)
H S ,H S , ∀k = 0, . . . , H S -1 .

Proof. See in Appendix E. Now we will provide an explicit upper bound for the operator norm L .

Proposition 9. The norm of the operator L = (L (1) , L (2) ), where

L (1) (X) = (L (1) 1 (X), . . . , L (1)
M (X)) , L (2) (X) = (L (2) 0 (X), . . . , L (2)
H S -1 (X)) , with L (1) m and L (2)
n2 defined in (41)-( 42), satisfies

L 2 L (1) 2 + L (2) 2 = T H S 2 + T M +1 2 = (H S + 1) + (M + 1) .
Proof. See in Appendix F.

To solve the problem [START_REF] Gelfand | Properties and operations[END_REF], we first propose Algorithm 1, which uses the primaldual method introduced in [START_REF] Polisano | Convex super-resolution detection of lines in images[END_REF]. Following [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]Theorem 5.1], we know that the method converges to a solution ( X, ξ0 , ..., ξQ-1 ) of the problem [START_REF] Hough | Method and means for recognizing complex patterns[END_REF], provided the parameters τ > 0 and σ > 0 in Algorithm 1 are such that

(55) 1 τ -σ L 2 > β 2 .
We then choose 0

< τ < 2, σ = (H S + M + 2) -2 (1/τ -β/1.9) and ρ n ≡ ρ = 1.
Algorithm 1 Primal-dual splitting algorithm for [START_REF] Hough | Method and means for recognizing complex patterns[END_REF] Input: y 1-D FFT of the blurred and noisy data image y Output: x solution of the optimization problem [START_REF] Condat | A 2-D spectral analysis method to estimate the modulation parameters in structured illumination microscopy[END_REF] under constraints [START_REF] Condat | Cadzow denoising upgraded: A new projection method for the recovery of Dirac pulses from noisy linear measurements[END_REF] 1: Initialize primal and dual variables to zero X 0 = 0, ξ i,0 = 0, ∀i ∈ 1, Q 2: for n = 1 to Number of iterations do 3:

X n+1 = prox τ G (X n -τ ∇F (X n ) -τ Q-1 i=0 L * i ξ i,n ), 4: X n+1 = ρ n Xn+1 -(1 -ρ n )X n , 5
:

for i = 0 to Q -1 do 6: ξ i,n+1 = prox σH * i (ξ i,n + σL i (2X n+1 -X n )), 7: ξ i,n+1 = ρ n ξi,n+1 + (1 -ρ n )ξ i,n , 8:
end for 9: end for We detail below the other terms in lines 3 and 6 of Algorithm 1, involving the computation of proximity operators and adjoint operators. For more details on convex analysis, monotone operator theory and proximal splitting methods, we refer the reader to [START_REF] Boyd | Convex optimization[END_REF][START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Parikh | Proximal algorithms[END_REF]. [START_REF] Fernandez-Granda | Super-resolution and compressed sensing[END_REF]. Then we have, for every m, n 2 :

Set x 0 = 1 H S H S -1 n2=0 x[0, n 2 ], set G = ι B , with B defined in
prox τ G ( x, q) =      x[0, n 2 ] = x 0 if x 0 c x[0, n 2 ] = c otherwise q[m, 0] = c if q[m, 0] > c
.

Let P C be the projection operator onto the cone of positive matrix C; by Moreau identity [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF],

prox σH * i (M) = M -σprox H i σ 1 σ M = M -P C (M) .
Finally, we need to compute, in line 3 of the algorithm, the adjoint operators L * i , where the operators L i are defined in ( 41)- [START_REF] Fatemi | Sampling and reconstruction of shapes with algebraic boundaries[END_REF]. The dual variables (ξ i,n ) i in Algorithm 1 refer to Hermitian Toeplitz matrices M

(1)

m ∈ T H S +1 or M (2) n2 ∈ T M +1 .
By definition, the adjoint operators give the images (z

(1) m , q (1) m ) = L (1) * m M (1)
m and (z

(2) n2 , q (2) n2 ) = L (2) * n2 M (2)
n2 . According to the definitions ( 41)-( 42), for a primal variable X = ( x, q) the operators L

(1)

m and L (1)
n2 act respectively on the (m + 1)-th row of the images ( x, q) and on the (n 2 + 1)-th column of the image x only; so we can easily see, concerning the adjoint operators, that q

(2) n2 = 0 and z (2) n2 (resp. z (1) m , q (1)
m ) is zero, except at the corresponding column index n 2 (resp. row index m), where z (2) n2 [:,

n 2 ] = T * M +1 M (2) n2 , (56)
(z (1) m [m, :],

q (1) m [m, :]) = T * H S M (1) m , (57)
with the expression of the adjoint operators T * M +1 and T * H S given in Proposition 8. Thus, the operations on X n = ( x n , q n ) before applying prox τ G consist in a gradient descent step X n -τ ∇F (X n ), followed by an update of all its rows and columns due to the terms -τ

Q-1 i=0 L * i ξ i,n
, whose expressions are provided by ( 56) and (57). 4.2. Second algorithm design. We can note that in Algorithm 1, τ must be smaller than 2/β, which is a limitation in terms of convergence speed. To overcome this issue, we subsequently developed a second algorithm, similar to Algorithm 1, but with the data fidelity term A xy Y activated through its proximity operator, instead of its gradient. We consider solving the optimization problem by an overrelaxed version [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF] of the Chambolle-Pock algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]:

(58) X = arg min X∈H {G(X) + H(L(X))} , with now G = 1 2 A • -y 2 Y which is proximable, Hx = Q i=0 H i x i with H i = ι C for i < Q, where L i = L (2) i when 0 i H S -1, L i = L (1) i-H S +1 when H S i Q -1 and H Q = ι B with L Q = Id. So now, L 2 H S + M + 3.
Let τ > 0 and σ > 0 such that τ σ L 2 = 1, then the primal-dual Algorithm 2, with F = 0 and weights ρ n ≡ ρ = 1.9, which is an over-relaxed version of the Chambolle-Pock algorithm, converges to a solution ( X, ξ0 , ..., ξQ-1 ) of the problem [START_REF] Hough | Method and means for recognizing complex patterns[END_REF] [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]Theorem 5.1].

The Algorithm 2 requires computing prox τ G . Since we have

p = prox τ G ( x) ⇔ x -p = ∇(τ G)(p) , (59)
⇔ x -p = τ A * (Ap -y) , ⇔ x + τ A * y = (I + τ A * A)p ,
the proximity operator has the following expression:

prox τ G ( x) = (I + τ A * A) -1 ( x + τ A * y) ,
Algorithm 2 Primal-dual splitting algorithm for [START_REF] Nacereddine | Robustness of Radon transform to white additive noise: general case study[END_REF] Input: y 1-D FFT of the blurred and noisy data image y Output: x solution of the optimization problem [START_REF] Condat | A 2-D spectral analysis method to estimate the modulation parameters in structured illumination microscopy[END_REF] under constraints [START_REF] Condat | Cadzow denoising upgraded: A new projection method for the recovery of Dirac pulses from noisy linear measurements[END_REF] 1: Initialize all primal and dual variables to zero 2: for n = 1 to Number of iterations do 3:

Xn+1 = prox τ G (X n -τ Q-1 i=0 L * i ξ i,n ),
4:

X n+1 = ρ n Xn+1 -(1 -ρ n )X n 5
:

for i = 0 to Q -1 do 6: ξi,n+1 = prox σH * i (ξ i,n + σL i (2X n+1 -X n )), 7
:

ξ i,n+1 = ρ n ξi,n+1 + (1 -ρ n )ξ i,n 8:
end for 9: end for for which we propose below two ways of computing the inverse.

We proved in Lemma 6 that A * y = G * y Ȟ and then

(I + τ A * A) x = x + P xQ, P = τ G * G, Q = ȞT Ȟ .
The square matrices P and Q are of size p = M + 1 and q = H S . We have to solve (I + τ A * A) x = z; that is, x + P xQ = z. This kind of system can be solved by the mean of the Kronecker Product (2) as:

x + P xQ = z ⇐⇒ (I pq + Q ⊗ P T )vec( x) = vec(z).
where vec( x) denotes the vectorization of the matrix x formed by stacking the columns of x into a single column vector and I pq + Q ⊗ P T is a matrix of size pq × pq which can be inverted, giving access to vec( x) and then to x. Finally, the operator prox τ G can be seen as a large matrix-vector product.

Another option consists in operating on the rows x m of x, since the operator A acts on them:

(I + τ A * A) x = z ⇐⇒ (I + | g m | 2 Q) x m = z m , ∀m = 0, . . . , M .
This time, the operator prox τ G involves performing M + 1 matrix-vector products of size q × q, which appears to be more efficient in practice.

Extended problem formulation.

We now consider a data image b containing lines with no angle restriction, which extends the previous case by relaxing the assumption made in Hypothesis 1. We can decompose this image into the sum of two images b = b 1 + b 2 , with b 1 (resp. b 2 ) containing vertical (resp. horizontal) lines, that is with angles in (-π/4, π/4] (resp. outside this range). We can also define x 1 of size (M + 1) × H S and x 2 of size W S × (P + 1) with W S = W + 2S and P = (H -1)/2 such as A x 1 = b 1 and à x 2 = b 2 , where g 2 = (0 P -S , h, 0 P -S ) and à denotes the operator which multiplies each column vector x 2 [:, n 2 ] by the corresponding Fourier coefficient g 2 [n 2 ] and convolves it with the filter h; that is,

à x 2 = ( G 2 x 2 ) * h with G 2 = diag( g 2 [0], . . . , g 2 [P ]
). We finally define the Hermitian symmetry operator S 1 (resp. S 2 ), which to each column v = (v 0 , v 1 , . . . , v M ) (resp. row v = (v 0 , v 1 , . . . , v P )) associates the symmetric extension (v * M , . . . , v 0 , . . . , v M ) (resp. (v * P , . . . , v 0 , . . . , v P )). Let X 1 = ( x 1 , q 1 ) and X 2 = ( x 2 , q 2 ) be the optimization variables, living in spaces

H 1 = X 1 × Q 1 and H 2 = X 2 × Q 2 . Let H = H 1 × H 2 , X = X 1 × X 2 and Q = Q 1 × Q 2 .
The data fidelity term is now:

F (X 1 , X 2 ) = 1 2 F -1 1 S 1 A x 1 + F -1 2 S 2 Ã x 2 -y 2 F = 1 2 A 1 x 1 + A 2 x 2 -y 2 F , with A 1 = F -1 1 S 1 A, A 2 = F -1 2 S 2 Ã, where F 1 (resp. F 2 )
is the Fourier transform with respect to the columns (resp. rows) and • F is the Frobenius norm.

Proposition 10. The gradient of F is (60) ∇F (X 1 , X 2 ) = 1 2 A 1 * (A 1 x 1 + A 2 x 2 -y) A 2 * (A 1 x 1 + A 2 x 2 -y) , which is Lipschitz-continuous of Lipschitz constant β = 1 min(W,H) . Proof. See in Appendix H.
The image x 1 keeps the same kind of constraints as in the Algorithm 1, which act similarly on the image x 2 in a rotated way; that is, we define

L (3) m (X 2 ) = T P +1 ( x 2 [m, :]) , (61)
L (4) n2

(X 2 ) = T W S (fliplr( x 2 [:, n 2 ]), fliplr(q 2 [:, n 2 ])) , (62)
where fliplr denotes a flip from left to right on each column of the matrix. The boundary constraints on x 1 and x 2 are respectively given by:

B 1 = ( x 1 , q 1 ) ∈ H 1 : x 1 [0, n 2 ] = x 1 [0, 0] c 1 , q 1 [m, 0] c 1 , (63)
B 2 = ( x 2 , q 2 ) ∈ H 2 : x 2 [m, 0] = x 2 [0, 0] c 2 , q 2 [P, n 2 ] c 2 . (64)
Likewise, the inner product on spaces X 2 and Q 2 are:

z 1 , z 2 X2 = W S -1 m=0 z 1 [m, 0]z 2 [m, 0] + 2Re P n2=1 W S -1 m=0 z 1 [m, n 2 ]z 2 [m, n 2 ] * , (65)
z 1 , z 2 Q2 = 2Re P n2=0 W S -1 m=0 z 1 [m, n 2 ]z 2 [m, n 2 ] * , (66)
and so the adjoint of the operators remain the same.

Like before, we define L

(3) (X 2 ) = (L (3) 0 (X 2 ), . . . , L (3)
W S -1 (X 2 )), also L (4) (X 2 ) = (L (4) 0 (X 2 ), . . . , L (4)
P (X 2 )) and L = (L (1) , L (2) , L (3) , L (4) ). It is easy to show that

L 2 L (1) 2 + L (2) 2 + L (3) 2 + L (4) 2 ,
(H S + 1) + (M + 1) + (P + 1) + (W S + 1) .

Finally, we have (67) ( X1 , X2 ) = arg min

(X1,X2)∈H 1 2 A 1 x 1 + A 2 x 2 -y 2 F + ι B1 (X 1 ) + M m=1 ι C (L (1) m (X 1 )) + H S -1 n2=0 ι C (L (2) n2 (X 1 )) + ι B2 (X 2 ) + W S -1 m=0 ι C (L (3) m (X 2 )) + P n2=1 ι C (L (4) n2 (X 2 )) .
4.4. Inpainting problems. We now consider the case in which a binary mask is applied on the data image, like in Figure 3. The corresponding linear operator, denoted by M, consists in elementwise multiplication of the matrix b with a binary matrix, whose zero coefficients are the indices of the pixels unavailable to observation. We have M * = M. The data fidelity term becomes F (X) = 1 2 MF -1 1 S 1 A xy 2 F , whose gradient can be expressed as previously, with β = 1/W (since F -1 1 = 1/W and M = 1). The constraints remain the same as in ( 28) and the method is also easily transposable to the extended setting of subsection 4.3.

At this point, the first part of the process has been completed; that is, the image x has been restored from the degraded image y. From this image, we can for instance reduce the blur by applying other filters g r and h r with smaller spread and visualize the resulting image b r , passing the solution x through this new blur operator A r ; that is, b r = A r x.

5.

Recovering the line parameters by a Prony-like method. In this section, we present the method that underlies the second step of this work (see Figure 3), namely the estimation of the line parameters, which is related to the spectral estimation field. We now focus on estimating the parameters (θ k , α k , η k ), which characterize the K lines, from the solution of the minimization problem x (27), extended by Hermitian symmetry to m = -M, . . . , -1 beforehand. This requires the use of a classical spectral estimation method [START_REF] Stoica | List of references on spectral line analysis[END_REF][START_REF] Stoica | Spectral Analysis of Signals[END_REF].

5.1. Sketch of the 1-D Prony TLS method. The recovering procedure hereafter, based on [START_REF] Rahman | Total least squares approach for frequency estimation using linear prediction[END_REF], is an extended method of the famous Prony method [START_REF] Prony | Essai experimental et analytique : sur les lois de la dilatabilité de fluides élastique et sur celles de la force expansive de la vapeur de l'alkool, à différentes températures[END_REF]. Let us sketch this method, which is based on an annihilating property [START_REF] Blu | The generalized annihilation property-A tool for solving finite rate of innovation problems[END_REF]. Let z = (z 0 , . . . , z N -1 ) be a complex vector, whose components are: [START_REF] Prony | Essai experimental et analytique : sur les lois de la dilatabilité de fluides élastique et sur celles de la force expansive de la vapeur de l'alkool, à différentes températures[END_REF]

z i = K k=1 d k e j2πf k i , ∀i = 0, . . . , N -1 , with d k ∈ C, f k ∈ [-1/2, 1/2
) the parameters to retrieve and N 2K + 1.

Let ζ k = e j2πf k , we introduce the annihilating polynomial

H(ζ) = K l=1 (ζ -ζ l ) = K l=0 h l ζ K-l with h 0 = 1.
Then, we can note that for all r = K, . . . , N -1:

(69) K l=0 h l z r-l = K l=0 h l K k=1 d k ζ r-l k = K k=1 d k ζ r-K k K l=0 h l ζ K-l k H(ζ k )=0 = 0 .
By rearranging the elements z i in a Toeplitz matrix P K (z) of size (N -K) × (K + 1) and rank K as follows:

(70)

P K (z) =    z K • • • z 0 . . . . . . . . . z N -1 • • • z N -K-1    , (69)
can be written with h = (h 0 , . . . , h K ) as:

(z * h)(r) = 0, ∀r = K, . . . , N -1 ⇐⇒ P K h = 0 .
Consequently, the method consists in finding a right singular vector h = (h 0 , . . . , h K ) of the matrix P K (z) associated to the singular value zero. From the SDV decomposition

P K (z) = V 1 ΣV * 2 , it corresponds to the (K + 1)-th column of V 1 , that is h = V 1 [:, K + 1]. Thus, the roots of the polynomial H(ζ) = K l=0 h l ζ K-l are the searched complex values ζ k = e j2πf k and then f k = arg(ζ k )/(2π).
The Algorithm 3 below described the procedure for estimating the frequencies:

Algorithm 3 Prony

Input: z = (z 0 , . . . , zN-1 ) ∈ C N a vector of form ( 68) possibly corrupted by noise, with N 2K + 1 and K the number of frequencies to retrieve. Output: { fk } K k=1 the estimated frequencies. Finally, the complex amplitudes can be retrieved as well, by writing (68) in matrix form z = Ud, where d = (d 1 , . . . , d K ) and the matrix U of size

N × K is (71) U = a(f 1 ) • • • a(f K ) =        1 • • • 1 e -j2πf1 • • • e -j2πf K e -j4πf1 • • • e -j4πf K . . . . . . . . . e -j2π(N -1)f1 • • • e -j2π(N -1)f K       
and we recover the amplitudes by least-squares approximation:

(72) d = (U * U) -1 U * z .
5.2. Procedure of the line parameters estimation. We start with the angle parameters estimation by applying a Prony-like method onto the rows tm of the output solution x, since from ( 23) and ( 24) one have theoritically [START_REF] Recht | Necessary and sufficient conditions for success of the nuclear norm heuristic for rank minimization[END_REF] t

m = x [m, :] = K k=1 d m,k a(f m,k ) T ,
whose frequencies to estimate are related to the angles by f m,k = tan θ k W m. The number of lines can be found either by a dedicated method like the approximate Prony method (APM) [START_REF] Potts | Parameter estimation for exponential sums by approximate Prony method[END_REF] or by evaluating the rank of the dual matrices T M +1 (x[:, n 2 ]) at the end of the algorithm, from the decreasing profile of their eigenvalues. The output solution of the algorithm x, whose stopping criteria is met when this one is sufficiently close to the exact solution xx 1, then remains some noise. In the noiseless case, the Prony method is always able to recover the frequencies with infinite precision, if the number of samples N is greater than 2K. But in our case, the estimate fm,k is affected by some uncertainty m,k ; that is, fm,k = f m,k + m,k , due to the instability of root finding in presence of noise. Then, we propose to estimate the angle parameters by applying K linear regression to the data { fm,k

} 1 m M since fm,k = tan θ k W m + m,k ,
which leads to an estimation of the slope tan θ k and then to the angle θ k . The errors m,k committed by evaluating the frequencies f m,k have an amplitude which depends on m. Indeed, for a small m, the frequencies f m,k = tan θ k W m, are close to each other and the Prony method fails to accurately determine the frequencies. Consequently, it is preferable to start the linear regression with the largest values m m 0 , in order to space the frequencies on the unit circle. We have to make sure that for large values of m, the two extremal frequencies, say f m,1 0 and f m,K 0, are not close to -π and π, respectively, at the same time, as this would violate the separation criteria. Then, the angles θk estimated from the K linear regressions are used to form the matrices Ũm and to obtain the complex estimated amplitudes dm , by solving least-squares linear systems [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF].

Remark 8. In the preliminary version of this work [START_REF] Polisano | Convex super-resolution detection of lines in images[END_REF], we proposed a simplistic method consisting in averaging the estimates

θm,k = arctan(W fm,k /m) = arctan(W f m,k /m + W m,k /m) ≈ θ m,k + W m,k /m ,
whose error is actually amplified by a factor W/m; this gives bad results, in particular for a small m. Consequently, the mean θk = 1 M M m=1 θm,k did not lead to a robust estimation of the angles θk .

The previous estimation process of the frequencies fm,k and the angle estimation are possible, as long as the sorting process of the frequencies fm,k is related to the corresponding angles

θ 1 • • • θ k • • • θ K for
all m, which allows us to perform the linear regression with respect to m. It would not be possible to do the same with fn2,k = (n 2 tan θk -ηk )/W to estimate the offsets η k , performing the Prony method on the columns, because the affine relation does not preserve the order (one can find n and n such that fn,k1 fn,k2 and fn ,k1 fn ,k2 ). The frequencies are not uniquely determined, as they belong to an interval of length greater than one fn2,k ∈ [-(H S + M )/W, (H S + M )/W ] and above all we would lose the correspondence between the fn2,k and the previous estimated angles θk , which compromises the estimation of the η k . Thus, the solution is to perform the Prony method on the estimated vectors ẽk = {ẽ m,k } M m m0 , noticing from (24) that

d m,k /|d m,k | = e -j2πη k /W m ,
which leads to the frequency estimation of a single sampled exponential. This way, the correspondence between the angles θk and the offsets ηk is preserved.

Finally, regarding the amplitudes c k , taking the modulus of the dm,k 's leads to inaccurate estimates, since they are computed from the solution x and their amplitudes have been shrunk, due to the choice of a parameter c < c to remove noise. Like often with variational methods, this yields a bias towards zero and yields images with a loss of contras; that is, the line intensities are globally smaller. That is why we propose to perform a refitting step, to recover the amplitudes in a better way. Given the estimated 2-D atoms 1 cos θk a 2D ( θk , ηk ), we evaluate the amplitudes αk by applying least-squares estimation to the noisy data y:

(74) (α 1 , . . . , αK ) = arg min α1,...,α K K k=1 α k Ãk -y 2 , Ãk = 1 cos θk Aa 2D ( θk , ηk ) .
The procedure is summarized Algorithm 4 and is illustrated in Figure 5:

11:

Compute the vector dm = ( dm,1 , . . . , dm,K ) T = ( Ũ * m Ũm ) -1 Ũ * m tm .

12:

Compute the normalized values ẽm,k = dm,k /| dm,k |. 13: end for 14: for k = 1 to K do 15:

Form the vector ẽk = (ẽ m,k ) M m m0 .

16:

Compute ηk ← W * Prony(ẽ k , 1). 17: end for 18: Evaluate the amplitudes {α k } K k=1 by solving the least-squares problem (74).

Other related works and further comments. Below we discuss other approaches to estimate line parameters and discuss their effectiveness:

• There exist sophisticated methods, called debiasing methods [START_REF] Brinkmann | Bias reduction in variational regularization[END_REF][START_REF] Deledalle | On debiasing restoration algorithms: applications to total-variation and nonlocal-means[END_REF][START_REF] Deledalle | Clear: Covariant least-square refitting with applications to image restoration[END_REF], which could be considered to recover the line intensities α k , instead of the least-squares method (74) used here.

• For recovering the offsets η k , we also could apply the Prony method on the middle line (n 2 = 0) of the image x, since from [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] we have [START_REF] Shukla | Sampling schemes for multidimensional signals with finite rate of innovation[END_REF] l and since the arguments of these exponentials are uniquely determined. Provided that the latter frequencies are well separated, it would give a better estimation of the offsets ηk , but the correspondence with the angles θk would be lost, as mentioned previously. To reconnect them, one could rely on the work of [START_REF] Potts | Parameter estimation for exponential sums by approximate Prony method[END_REF] (or equivalent [START_REF] Plonka | How many Fourier samples are needed for real function reconstruction?[END_REF]), noticing that

0 [m] = K k=1 c k e j2πη k m/W
K k=1 d 0,k a(f 0,k ) . . . K k=1 d m,k a(f m,k ) . . . K k=1 d M,k a(f M,k ) t 0 t m t M x f0,k . . .
x [m, n 2 ] = K k=1 α k cos θ k e j2π(n2 tan θ k -η k ) m W = K k=1 α k cos θ k e j2πm ( n 2 -1 )( tan θ k /W η k /W )
are the samples (indexed by m) of the Fourier transform of

K k=1 α k cos θ k δ (tan θ k /W,η k /W ) (denoted by h 0 ) along radial lines indexed by n 2 , that is h 0 (mn 2 , m) = x [m, n 2 ].
One could exploit this information to couple the right frequencies (tan θ k /W, η k /W ) k among all the possibilities {(tan θ k1 /W, η k2 /W ) : k 1 , k 2 = 1, . . . , K}, using the estimated frequencies along the radial lines h 0 (mn 2 , m); that is, using the Prony method along columns of x.

This procedure, called sparse approximate Prony method (SAPM), requires the use of samples along lines h 0 (n, 0), h 0 (0, n), h 0 (n, αn + β) for n ∈ Z N = [-N, N ] ∩ Z and α, β ∈ Z, conveniently chosen. Therefore, it first requires to separately estimate the frequencies along the horizontal axis h 0 (n, 0) and the vertical axis h 0 (0, n), before coupling them through the lines h 0 (n, αn + β). We get the frequencies η k /W along one axis, since from [START_REF] Shukla | Sampling schemes for multidimensional signals with finite rate of innovation[END_REF] we have h 0 (0, m) = l 0 [m] = x [m, 0]. However, we do not directly have access to h 0 (n, 0) from x for estimating the tan θ k /W . We admittedly have these frequencies appearing on the first row x [1, :], but they are clustered and badly separated on this one. Consequently, it is better to use our procedure for estimating the angles, which exploits the other rows for larger index m, in order to space the frequencies before applying the Prony method; this is reminiscent of the strategies of decimation developed in [START_REF] Cuyt | Multivariate exponential analysis from the minimal number of samples[END_REF][START_REF] Briani | VEXPA: Validated EXPonential Analysis through regular sub-sampling[END_REF][START_REF] Batenkov | Accurate solution of near-colliding Prony systems via decimation and homotopy continuation[END_REF]. Moreover, our proposed method to estimate the offsets η k deals with the rows and columns jointly and automatically preserves the correspondence, without having to restore it a posteriori. This is less arbitrary than performing the Prony method on l 0 (75), which requires well-separated frequencies on this column; this is a meaningless condition, since the lines can intersect the x-axis.

• The authors of [START_REF] Dossal | Sampling the Fourier transform along radial lines[END_REF] proposed a convex approach to the recovery of a superposition of point sources from samples of its Fourier transform along radial lines. They also emphasize the equivalence between working with the radial Fourier coefficients and working with its Radon projections, according to the Fourier slice theorem. An important question is then how many lines (or projections) and samples are needed to obtain an exact reconstruction. The authors of [START_REF] Renyi | On projections of probability distributions[END_REF] show that K points can theoretically be determined by the projection onto K + 1 distinct lines through the origin. In [START_REF] Plonka | How many Fourier samples are needed for real function reconstruction?[END_REF], it was conjectured that under certain additional assumptions, it is possible to choose only four lines passing through the origin to guarantee a unique reconstruction, but this conjecture has been shown to be false in [START_REF] Diederichs | Parameter estimation for bivariate exponential sums[END_REF]. These authors cast the problem into a nonconvex optimization problem for which there are guarantees for recovering the points, but this constrained optimization problem is NP-hard. Finally, the authors of [START_REF] Dossal | Sampling the Fourier transform along radial lines[END_REF] give a total variation minimization formulation of the problem. Like in [START_REF] Plonka | How many Fourier samples are needed for real function reconstruction?[END_REF], where the radial line of angle θ is chosen to maximize the minimal distance between two projections onto it, the authors show the crucial role of the minimal separation distance ν min (∆ θ ) between the projected positions ∆ θ = { θ, x k } K k=1 . Thus, they proved that for a set of radial measurements Θ ∈ S 1 for which the global minimal separation distance is ν min = inf θ∈Θ ν min (∆ θ ), total variation minimization has a unique solution achieved by their algorithm from 3 different radial lines of Θ and a sampling set Z N with N 2/ν min . These results provide strong guarantees, but the framework is quite different from ours, especially because the unit sampling step is constant on every radial line, while in our case the samples are taken on points (mn 2 , m). Moreover, the frequencies along these lines can easily collide, which contradicts the separation assumption on the projected positions. Finally, a common point between our two-step super-resolution recovery of lines and the references cited above is that all these methods come down to 1-D settings, based on a sufficient separation property.

• A final approach might be to apply the Radon transform to the noisy blurred lines y = b + , which somehow transform blurred lines to blurred peaks, representing in Figure 1b by 2-D point sources in the Radon domain. Then, we come down to a classical 2-D spikes super-resolution problem, with on the one hand a noise amount which is somehow reduced [START_REF] Nacereddine | Robustness of Radon transform to white additive noise: general case study[END_REF] and on the other hand a minimal separation inf k =l dist(θ k , θ l ), which is better than the one we have with the frequencies tan θ k /W of x. However, the main drawback is that the Radon transform is performed on the data image y of blurred and noisy lines, which is finite and discrete, so that it produces some artifacts like "Butterfly patterns" (see Figure 1) [START_REF] Davies | Simple two-stage method for the accurate location of Hough transform peaks[END_REF][START_REF] Lassen | Automated Determanation of Crystal Orientations from Electron Backscattering Patterns[END_REF][START_REF] Leavers | The Radon transform and its application to shape parametrization in machine vision[END_REF], due to the finite length of lines. Moreover, they induce discrete other approximations due to the projections on a grid [START_REF] Vasconcelos | On the Radon transform of sampled functions[END_REF]. These problems are avoided with our model and assumptions, which yields an off-the-grid super-resolution estimation of the parameters and outperforms the naive approaches using the Radon (or Hough) transform.

7. Experimental results. The reconstruction procedure was implemented in matlab™ code, available on the webpage of the first author. We consider an image of size W = H = 65, containing three lines of parameters (θ 1 , η 1 , α 1 ) = (-π/5, 0, 255), (θ 2 , η 2 , α 2 ) = (π/16, -15, 255) and (θ 3 , η 3 , α 3 ) = (π/6, 10, 255). We consider the normalized filter h approximating a Gaussian function of standard deviation κ; that is, ϕ 2 : t → (2πκ 2 ) -1/2 exp(-t 2 /(2κ 2 )), on the compact set [-S, S] with S = 4κ -1; Relative errors Fig. 6. Log-log plot of the relative errors

x-x X x X and A x-y Y y Y
for the first experiment.

and the normalized filter g = (0 M -S , h, 0 M -S ), whose DFT is an interpolation of h, which approaches the continuous Fourier transform

ϕ 2 : ν → exp(-2π 2 κ 2 ν 2 ). Then, g ∞ = h ∞ = 1.
We use the Algorithm 2 for solving the optimization problem [START_REF] Gelfand | Properties and operations[END_REF] and recovering the lines in the following experiments:

• The first experiment consists in the reconstruction of the lines from x in absence of noise, (1) by applying the operator A on this solution, possibly with other kernels g r and h r and then taking the 1-D inverse Fourier transform; and (2) by applying the Prony method to recover the parameters of the lines, in the aim to display the lines by vectorial drawing. We run the algorithm for 10 6 iterations. Results of relative errors for the solution x and the estimated parameters are given in Figure 6 and Table 1, where

∆ θi /θ i = |θ i -θi |/|θ i |, ∆ αi /α i = |α i -αi |/|α i | and ∆ ηi = |η i -ηi |.
Although the algorithm is slow to achieve high accuracy, convergence is guaranteed and we observe empirically perfect reconstruction of x , when the lines are not too close to each other.

• The purpose of the second experiment is to show the robustness of the method in presence of a strong noise level (Figure 7a). With c = c /3 and only 2.10 3 iterations, we are able to completely remove the noise and to estimate the line parameters with an error of 10 -2 .

• Finally, the third experiment, with 10 5 iterations, illustrates the efficiency of the method even in presence of a large blur (Figure 7b), yielding an error of 10 -4 . For both experiments, the estimated images corresponding to step (1) and ( 2) are visually identical and are displayed in Figure 7c andFigure 7d.

We emphasize that our algorithm has a much higher accuracy than what is achieved by peak detection after the Hough or Radon transform. These methods are relevant for giving a coarse estimation of line parameters. They are robust to strong noise, but completely fail with a strong blur, which prevents peak detection (see Figure 8). Notice that even by decreasing the discretization steps of the process, we rapidly reach a plateau, as illustrated by Figure 9. This method is limited in accuracy by the pixel grid. By contrast, our super-resolution method makes it possible to achieve infinite precision for the line parameters. ∆ η (4.10 -6 , 7.10 -6 , 7.10 -6 ) (5.10 -2 , 4.10 -2 , 3.10 -2 ) (5.10 -5 , 10 -4 , 3.10 -4 )

7.1. Closing lines. For a reasonable amount of noise (ζ = 20), the algorithm succeeds in separating two close lines (θ 1 , η 1 , α 1 ) = (-0.73, -1, 255) and (θ 2 , η 2 , α 2 ) = (-0.75, 1, 255) as illustrated in Figure 10. The estimation of the parameters gives ( θ1 , η1 , α1 ) = (-0.725, -0.7, 237) and ( θ2 , η2 , α2 ) = (-0.753, -0.6, 251).

More lines and different amplitudes.

A more complicated example is depicted in Figure 11a, containing seven well-separated lines whose parameters are enumerated in Table 2, corrupted by some noise with variance ζ = 20. We run the algorithm with c = 0.8c , τ = 1, σ = (τ (M +H S +2)) -1 and after only 2.10 3 iterations, we are able to denoise the image as illustrated in Figure 11b and to estimate the line parameters, with the proposed Prony procedure, as illustrated in Figure 11d, with an error of 10 -2 , as reported in Table 3.

General case.

We consider an image y = b + composed as four noisy blurred lines (κ = 1 and ζ = 20) whose two of them are rather vertical (i.e with angle in (-π/4, π/4]), while the two other ones are rather horizontal. The extended algorithm presented in subsection 4.3, provides after n = 6 • 10 4 iterations the denoised images depicted in Figure 12. It acts as an angles selector according to the horizontality and verticality: the rather vertical lines are gathered in the reconstructed image b1 (see Figure 12b) and the rather horizontal lines are gathered in the reconstructed image b2 (see Figure 12c), so that b ≈ b1 + b1 , which validates the procedure. 7.4. Inpainting. We consider occluded lines with a mask M as described in the subsection 4.4. In Figure 13, we occlude a portion of the blurred line (κ = 1) by applying a square mask M in the middle. We run the optimization algorithm and visualize the Figure 13a to Figure 13d, that is the evolution of the reconstruction after n = 2.10 3 , n = 10 4 and n = 10 6 iterations. We can see that the method is able to reconstruct the part of line occluded thanks to the information available outside the mask and to the optimization constraints related to the line structure. In the same way, Figure 13e to Figure 13h enable to visualize the "diffusion process" of the information within the occluded part in the Fourier domain. Finally, the experiment is carried out for an image containing three occluded lines, either by a bigger mask in Figure 14a, or by a mask whose binary inputs were drawn randomly Figure 14b. In both cases the lines are successfully reconstructed as in Figure 14c.

Conclusion.

We proposed a new variational formulation for the problem of recovering lines in degraded images, using the framework of atomic norm minimization. A primal-dual splitting algorithm has been used to solve the convex optimization problem. We applied it successfully to several image restoration problems, recovering lines parameters by the Prony method and we showed the robustness of the method to strong blur and strong noise level. We insist on the novelty of our approach, which is to estimate lines with parameters (angle, offset, amplitude) living in a continuum, with perfect reconstruction in absence of noise, without being limited by the discrete nature of the image, nor its finite size. This work can be viewed as a proof of concept for super-resolution line detection and invites us to revisit the Hough transform in a continuous way. Many theoretical questions remain open, like the study of the separation conditions under which perfect reconstruction can be guaranteed. The robustness of the method needs to be theoretically studied and would requires a statistical analysis ; this is left for future work. From a practical point of view, parallel computing would be welcome to speed up the proposed algorithm. We should also investigate the possibility of relaxing the periodicity and bandlimitedness assumptions, possibly by solving a convex feasibility program without any regularizers [START_REF] Eftekhari | Sparse non-negative super-resolution-simplified and stabilised[END_REF][START_REF] Eftekhari | Stable super-resolution of images[END_REF]. At the time of finalizing this paper, we became aware of the two other recent papers [START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF] and [START_REF] Eftekhari | A bridge between past and present: Exchange and conditional gradient methods are equivalent[END_REF]; we leave for future work the study of their relationship to the setting considered here. In subsequent work, we plan to apply the proposed approach to biomedical images containing curved structures, like tubulins, by operating on small overlapping patches. (78)

ϕ 1 (t) = 1 W W -1 n=0 g[n]D M 2π(t -nW ) W , ∀t ∈ R .
Proof. Let us start with the following classical theorem; in our notations:

Theorem 12 (Nyquist Whittaker-Shannon). Let f be a T -periodic function and c m (f ) = 0 for |m| M + 1. Then f can be reconstructed from the regular sampling {f (ka), k = 0, 1, ..., 2M }, where a = T 2M +1 is the sampling rate, in this way:

(79) f (x) = 1 2M + 1 2M k=0 f (ka)D M 2π T (x -ka) ,
where D M is the Dirichlet Kernel:

(80) D M (x) = M m=-M e jmx = sin (M + 1 2 )x sin x 2 .
Applying to ϕ 1 which is W -periodic and c m (ϕ 1 ) = 0 for |m| M +1 with 2M +1 = W , then a = 1 and the Nyquist theorem leads to:

(81) ϕ 1 (t) = W -1 n=0 ϕ 1 (n) sin π(t -n) W sin π(t-n) W .
Consequently, we give this explicit formula g[n] = ϕ 1 (n) and g[n] = c n-M (ϕ 1 ) for n = 0, 1, ..., 2M . By the simple change of variable n ← n -M , one can also obtain:

(

) M n=-M ϕ 1 (n) sin π(t -n) W sin π(t-n) W . Then, 1 W W 0 ϕ 1 (t) dt = 1 W M n=-M g[n] = 1 , 82
that is, the filter g is normalized as well. Moreover

W 0 ϕ 1 (t) 2 dt = M n=-M g[n] 2 = M n=-M ϕ 1 (n) 2 ,
and by Parseval relation:

m∈Z |c m (ϕ 1 )| 2 = M m=-M | g[m]| 2 = 1 W W 0 ϕ 1 (t) 2 dt .
Now, we describe the sampling process leading from continuous to discrete formulation, based on the following proposition:

Proposition 13. It is equivalent to perform the vertical convolution of u = x * φ 1 with ϕ 2 , with ϕ 2 * sinc, or with the Dirac comb γ :

t 2 → S n=-S h[n]δ(t 2 -n), where h[n] = (ϕ 2 * sinc)(n).
Proof. Due to Hypothesis 1, the assumption θ k ∈ (-π/4, π/4] yields | tan θ k | 1 for every k = 1, . . . , K. So, the function u given in [START_REF] So | Linear prediction approach for efficient frequency estimation of multiple real sinusoids: algorithms and analyses[END_REF], as a function of t 2 at fixed t 1 , is bandlimited: for every t 1 ∈ [0, W ), the Fourier transform F 2 u : ω 2 → R u (t 1 , t 2 )e -j2πω2t2 dt 2 , which is a distribution (sum of K Dirac combs), is zero for every |ω 2 | 1/2. Indeed, we have:

[F 2 u ](ω 2 ) = K k=1 α k sin θ k ϕ 1 ω 2 tan θ k exp j2πω 2 t 1 -η k tan θ k . Since | tan θ k | 1, we have |ω 2 / tan θ k | |ω 2 |. The support of ϕ 1 is included in [-1/2, 1/2]
(c m (ϕ 1 ) = 0 for |m| M + 1 and M/W < 1/2), as well as the support of F 2 u which is necessarily included in the support of ϕ 1 . Then, we have the equivalence

F 2 u = F 2 u • 1 [-1/2,1/2] ⇔ u = u * sinc and furthermore u * ϕ 2 = u * (ϕ 2 * sinc). In the Fourier domain, the function h = ϕ 2 * sinc is bandlimited, so [F 2 u ] h = [F 2 u ] h per
where h per corresponds to the periodization of the spectrum of h with period 1, which amounts to saying that

u * h = u * n h[n]δ(• -n) .
Remark 10. Assumption (ii') implies that the filter (h[n]) n should have compact support, but we can note that the function h = ϕ 2 * sinc does not have compact support, since it is bandlimited. This means that the continuous function h has to vanish at integer points t = n for |n| > S. Given such a compact filter (h[n]) S n=-S , the unique bandlimited function h satisfying these conditions is obtained by the Shannon interpolation formula:

h(t) = S n=-S h[n] sinc(t -n) .
By uniqueness, we necessarily have ϕ 2 * sinc = h and we can notice that there always exists a bandlimited solution ϕ 2 of this equation, which is simply ϕ 2 = h. In practice, we can always approximate a PSF by a bandlimited function h, with 2S + 1 samples h[n] of this PSF; that is why we argued in Assumption (ii) that the compact support assumption is not restrictive. Now, to obtain the discrete image b of ( 5), let us first define u by sampling u = x * φ 1 with unit step:

(83) u [n 1 , n 2 ] = (x * φ 1 )(n 1 , n 2 ), ∀n 1 = 0, . . . , W -1, n 2 = -S, . . . , H -1 + S .
With the above assumptions and Proposition 13, we can express b from u using a discrete vertical convolution with the filter h:

(84) b [n 1 , n 2 ] = S p=-S u [n 1 , n 2 -p]h[p], ∀n 1 = 0, . . . , W -1, n 2 = 0, . . . , H -1 . T N (a(f k )) =      1 e -j2πf k • • • e -j2πf k (N -1) e j2πf k 1 • • • e -j2πf k (N -2) . . . . . . . . . . . . e j2πf k (N -1) e j2πf k (N -2) • • • 1      , =      1 e j2πf k . . . e j2πf k (N -1)      1 e -j2πf k • • • e -j2πf k (N -1) , = a(f k )a(f k ) * .
We deduce that

T N (q) = K k=1 c k T(a(f k )) , = K k=1 c k a(f k )a(f k ) * , = K k=1 c k a(f k , φ k )a(f k , φ k ) * .
Therefore, the matrix

T N (q) z z * q 0 = K k=1 c k a(f k , φ k ) 1 a(f k , φ k ) 1 * is positive semidefinite. Given q 0 = K k=1 c k , we get SDP(z) K k=1 c k .
Since this holds for any decomposition of z, we conclude that SDP(z) z A .

• Conversely, let q ∈ C N be a vector such that q 0 0 and T N (q) z z * q 0 0.

In particular we have T N (q) 0. We denote by r the rank of T N (q). Theorem 2 ensures that T N (q) 0 and is of rank r N , if and only if there exists d k > 0 and distinct f k , such that

q = r k=1 d k a(f k ) , (89)
q 0 = r k=1 d k . ( 90
)
Let us set D = diag(d 1 , . . . , d r ) and

V = a(f 1 ) • • • a(f r ) =        1 1 • • • 1 e j2πf1 e j2πf2 • • • e j2πfr e j2πf12 e j2πf22 • • • e j2πfr2 . . . . . . . . . . . . e j2πf1(N -1) e j2πf2(N -1) • • • e j2πfr(N -1)       
.

By linearity of the operator T N :

T N (q) = r k=1 d k T N (a(f k )) , = r k=1 d k a(f k )a(f k ) * , = VDV * .
Since T N (a(f k )) contains only ones on the diagonal, we have

1 N tr(T N (q)) = r k=1 d k > 0 .
Besides, 1 N tr(T N (q)) = q 0 , therefore q 0 > 0. Let be M a general block matrix M = A B B * C , the Schur complement gives

[ C 0 ⇒ M 0 ] =⇒ [ A -BC -1 B * 0 ] .
We apply this lemma to M = T N (q) z z * q 0 , with A = T N (q), B = z and C = q 0 .

The left term is satisfied by hypothesis, hence

T N (q) -q -1 0 zz * 0 ⇐⇒ VDV * -q -1 0 zz * 0 . We define the square matrix V r by extracting the r first rows and columns of V, which is a Vandermonde matrix, whose determinant is det(V r ) = 1 k<l r (a(f l ) -a(f k )) .

Since we assumed f k = f l , ∀k = l, V r is invertible and rank(V) = r. Let us define v : C r → C N and v * : C N → C r the linear operators corresponding to matrices V and V * . We have rank(v * ) = rank(v) = r. By the rank-nullity theorem: dim(ker v * ) = N -r .

Thus, there exists a vector p ∈ C N such that p = 0 and V * p = 0 ⇔ p * V = 0. Consequently, p * (VDV * -q -1 0 zz * )p 0 ⇔ (p * V)D(V * p) -q -1 0 p * zz * p 0 , ⇔ q -1 0 p * z 2 by definition of the atomic norm [START_REF] Chandrasekaran | The convex geometry of linear inverse problems[END_REF]. Taking the infimum leads to SDP(z) z A .

• Finally, let us show that the infimum of the linear form : q → q 0 is achieved on the set A(z) = q ∈ R + × C N -1 : T N (z, q) = T N (q) z z * q 0 0 ; that is, (91) SDP(z) = inf q∈A(z)

(q) = min q∈A(z) (q) .

(i) Let us notice that since q ∈ A(z) implies T N (q) 0 then (q) = q 0 = 1 N tr(T N (q)) = N -1 i=0 λ i 0 , with λ i the eigenvalues of T N (q) which are positive reals.

(ii) First we show that A(z) is nonempty, since q = ( z 2 , 0, . . . , 0) T ∈ A(z). Indeed for a fixed vector z = (z 0 , . . . , z N -1 ) ∈ C N , v = (v 0 , . . . , v N ) ∈ C N +1 and v = (v 0 , . . . , v N -1 ) ∈ C N we have for this q :

v * T N (q) z z * q 0 v = z 2 v 2 2 + 2 Re v N N -1 i=0 z i v * i , z 2 v 2 2 -2|v N || z, v | , z 2 v 2 2 -2|v N | z 2 v 2 , z 2 ( v 2 2 -2|v N | v 2 + |v N | 2 ) , z 2 ( v 2 -|v N |) 2 , 0 .
Then q = ( z 2 , 0, . . . , 0) T ∈ A(z) and q 0 = z 2 , which means that A(z) is nonempty and the set { (q) : q ∈ A(z)} ⊂ R + is nonempty, so it admits a lower bound less than or equal to z 2 , hence 0 SDP(z) = inf q∈A(z) (q) z 2 .

(iii) From (ii), we have

SDP(z) = inf q∈A(z) (q) = inf B(z) (q) ,
where B(z) = {q ∈ A(z), q 0 z 2 } ⊂ A(z) . Now, from (89) and (90) we can show that B(z) is bounded since:

∀q ∈ A(z), q 2 r k=1 d k a(f k ) 2 √ N r k=1 d k = √ N q 0 , hence ∀q ∈ B(z), q 2 √ N q 0 √ N z 2 .
Consequently, B(z) ⊂ B • 2 (0, z 2 ) .

(iv) Moreover, A(z) = T N (z, •) -1 (C) is a closed set since the cone of positive matrix C is closed and the mapping T N (z, •) is linear, so it is continuous in finite dimension. Thus, B(z) = A(z) ∩ q ∈ R + × C N -1 : q 0 z 2 is a closed set, as the intersection of two closed sets.

(v) From (iii) and (iv), we conclude that B(z) is a compact set and that the function , which is linear and then continuous, achieves its minimum on B(z) so on A(z), which proves the result (91).

Appendix D. Proof of Proposition 4.

The proof of the direct implication is straightforward. Let us consider the converse one.

By Theorem 2, since ∀n, T M (l n ) 0 and is of rank one, then there exists γ n 0 and f n ∈ [0, 1[ such that Then, by writing

T N (z, q) =      0 T N (q) . . . 0 0 • • • 0 q 0      +      z 0 0 . . . z N -1 z * 0 • • • z * N -1 0     
we obtain as well for M ∈ T N +1 :

T N (z, q) , M M = 2Re Consequently, the adjoint of the operator

T M +1 : (X l , • , • X l ) → (T M +1 , • , • M ) ,
with the inner product • , • X l defined in [START_REF] Donoho | Compressed sensing[END_REF], when applied to a matrix M (2) ∈ T M +1 , yields the vector z 2 = T * M +1 M (2) ∈ R × C M , whose components are:

z 2 [k] = M -k l=0 M (2)
l+k,l , ∀k = 0, . . . , M .

Similarly, the adjoint of the operator

T H S : (X t × Q t , • , • Xt + • , • Qt ) → (T H S +1 , • , • M ) ,
with the inner products • , • Xt and • , • Qt defined in (36)- [START_REF] Dragotti | Sampling moments and reconstructing signals of finite rate of innovation: Shannon meets Strang-Fix[END_REF], when applied to M (2) ∈ T H S +1 , yields the pair of vectors (z 1 , q 1 ) = T * H S M (1) ∈ C H S × (R × C H S -1 ) , that is (dF 1 ) x (h) = A * (A xy), h X .

Moreover, F 2 is linear so (dF 2 ) X0 (X) = F 2 (X), hence (dF ) X0 (X) = A * (A x 0 -y), x X = A * (A x 0 -y) 0 , X H , that is ∇F (X 0 ) = A * (A x 0 -y) 0 .

Consequently,

∇F (X) -∇F (X ) H A * A( x -x ) X A * A x -x X .
We get β = A * A = A 2 and Lemma 5 concludes the proof.

We have A = A * = 1, S i x 1 F = x i X , for i ∈ {1, 2}; that is, S i = 1 and

F -1 i v 2 2 = 1 N 2 v 2 2 ; that is, F -1 i = 1 N . Hence, A 1 1 W , A 1 * 1, A 2 1 H
and A 2 1. Consequently, we get

C 1 1 2W 2 x 1 -x 1 2 X1 + 1 2H 2 x 2 -x 2 2 X2
and exactly the same majoration for C 2 . Thus, we have

∇F (X 1 , X 2 ) -∇F (X 1 , X 2 ) 2 X β 2 ( X 1 -X 1 2 H1 + X 2 -X 2 2 H2 ) , with β = 1 min(W, H)
.

Fig. 1 .

 1 Fig. 1. (a) The image b of three blurred lines and (b) the Radon transform of b .

Fig. 2 .

 2 Fig. 2. (a) Parameters (θ, η) characterizing the position of a line in the 2-D plane, (b) the matrix convention we use to display the image obtained by sampling with unit step the blurred line b = x * φ and (c) the resulting discrete image b [n 1 , n 2 ] = (x * φ)(n 1 , n 2 ).

2. 3 .

 3 Toward an inverse problem in the Fourier domain. Let us further characterize the blurred image b in Fourier domain, i.e. its 1-D discrete Fourier transform (DFT) on its first component (in the horizontal direction) denoted by b .

Fig. 3 .

 3 Fig.3. The two steps of the procedure : a convex optimization formulation for the reconstruction of the lines (in orange) and a Prony-like method for the estimation of their parameters (in purple).

1 : 2 : 3 :

 123 Compute the SVD decomposition Ṽ1 Σ Ṽ * 2 of the matrix P K (z) (70) Extract the (K + 1)-th right singular vector h = ( h0 , . . . , hK ) T = Ṽ1 [:, K + 1]. Compute fk = arg( ζk )/(2π), with { ζk } K k=1 the roots of polynomial K k=0 hk ζ k .

Algorithm 4 3 :

 43 Extraction of line parameters Input: x solution of the optimization problem (27) under constraints (28) Output: { θk } K k=1 , {η k } K k=1 and {α k } K k=1 the estimated parameters of the lines x . 1: Estimate K via APM or the rank of the matrices T M +1 (x[:, n 2 ]). 2: for m = 1 to M do Extract the row tm = x[m, :].

4 : 7 : 10 :

 4710 Compute { fm,k } K k=1 ← Prony( tm , K) and sort the frequencies with resp. to k. 5: end for 6: for k = 1 to K do Perform a linear regression on { fm,k } M m m0 to estimate tan θk and then θk . 8: end for 9: for m = 1 to M do Form the matrix Ũm = (a(tan θ1 m/W ) • • • a(tan θK m/W )).

Fig. 5 .

 5 Fig. 5. Parameters extraction procedure illustrating the Algorithm 4.

Fig. 7 .

 7 Fig. 7. (a) An image y = b + of three lines x , blurred by a Gaussian kernel φ with spread κ = 1 and corrupted by a strong noise ∼ N (0, ζ) with ζ = 200, for the second experiment, (b) with a strong blur (κ = 8) and no noise (ζ = 0) for the third experiment, (c) the denoised image b, (d) a vectorial drawing of the estimated lines of x by the Prony-like method.

Fig. 8 .

 8 Fig. 8. The Radon transform of the image y for Experiments 1, 2 and 3. The true parameters of the lines are in green.

16 Fig. 9 .Fig. 10 .

 16910 Fig. 9. Change of accuracy with respect to the angle step of the matlab™ Radon transform.

Fig. 11 .

 11 Fig. 11. (a) An image y = b + of seven well-separated lines x , blurred by a Gaussian kernel φ with spread κ = 1 and corrupted by noise ∼ N (0, 20) (b) the denoised image b, (c) the ground truth image b (d) a vectorial drawing of the estimated lines of x by the Prony-like method.

∆θ 1 .

 1 10 -2 2.10 -2 1.10 -3 2.10 -3 5.10 -3 5.10 -3 1.10 -3 ∆η 5.10 -1 7.10 -2 4.10 -2 1.10 -1 1.10 -2 2.10 -2 1.10 -2 ∆α/α 4.10 -2 5.10 -2 5.10 -3 4.10 -2 6.10 -3 1.10 -2 4.10 -3

Fig. 12 .Fig. 13 .

 1213 Fig. 12. (a) An image y = b + of four well-separated lines x , blurred by a Gaussian kernel φ with spread κ = 1 and corrupted by noise ∼ N (0, 20), two of them being rather vertical (angle in (-π/4, π/4]) and the two other ones rather horizontal, (b) the reconstructed image b1 , (c) the reconstructed image b2 .

Fig. 14 .

 14 Fig. 14. (a) An image y = b + of three line x , blurred by a Gaussian kernel φ with spread κ = 1, no noise and a large square mask M, (b) with a random mask and (c) their reconstruction at convergence.

Fig. 15 .

 15 Fig. 15. Relations between the variables.

  ⇔ p * z = 0 , ⇔ p ⊥ z .Since p ∈ ker v * , then z ∈ (ker v * ) ⊥ = Im v, so there exists a vector w ∈ C r such that z = Vw = r k=1 w k a(f k ), henceVDV * -q -1 0 Vww * V * 0 . Besides, Im v * ⊂ C r and dim(Im v * ) = rank(v * ) = r = dim(C r ), thus Im v * = C r and v * is surjective. Consequently, there exists a vector u ∈ C N such that V * u = sgn(w) = (w 1 /|w 1 |, . . . , w r /|w r |) T and u * (VDV * -q -1 0 Vww * V * )u 0 ⇔ (u * V)D(V * u) -q -1 0 (u * V)ww * (V * u) 0 ,⇔ sgn(w) * Dsgn(w) -1 q 0 sgn(w) * ww * sgn(w) 0,

P 1 (M

 1 l n [m] = γ n exp(j2πf n m) .Since we assume ∀n, x[0, n] = x[0, 0] = c 1 , then l n [0] = l 0 [0] = c 1 , we have (92) l n [m] = c 1 exp(j2πf n m) .Let m be fixed. The Prony matrix P 1 (t m ) of size 2 × (N -1)t m [N -1] t m [N -2]   is of rank one, consequently there exists λ m ∈ C such thatt m [n + 1] = λ m t m [n], ∀ 0 n N -2 .Thus,t m [n] = λ n m t m [0], ∀ 0 n N -1 . From (92) t m [0] = l 0 [m] = c 1 exp(j2πf 0 m), t m [1] = l 1 [m] = c 1 exp(j2πf 1 m) and then λ m = t m [1] t m [0] = 1 [m] 0 [m] = exp(j2π(f 1 -f 0 )m) .Therefore, we havet m [n] = λ n m t m [0] , = exp(j2π(f 1 -f 0 )m) n c 1 exp(j2πf 0 m) , = c 1 exp[j2π((f 1 -f 0 )n + f 0 )m] .Appendix E. Proof of Proposition 8. For z = (z 0 , . . . , z N -1 ) ∈ R × C N -1 and M ∈ T N ⊂ M N a Hermitian Toeplitz matrix of dimension N . We have:T N (z) , M M = 0 i,j N -1 [T N (z)] * ij M ij , = 0 i j N -1 z j-i M ij + 0 j<i N -1 z * i-j M ij , l+k,l,with ( * ) a change of variable k ← j -i and ( * * ) using that M l,l+k = M * l+k,l .

Table 1

 1 Errors on line parameters recovered by the proposed method.(10 -7 , 3.10 -6 , 7.10 -7 ) (10 -2 , 6.10 -2 , 9.10 -2 ) (6.10 -7 , 9.10 -5 , 8.10 -6 )∆ α /α (10 -7 , 10 -7 , 10 -7 ) (10 -2 , 9.10 -2 , 2.10 -1 ) (4.10 -5 , 2.10 -5 , 2.10 -5 )

	Experiment 1	Experiment 2	Experiment 3
	∆ θ /θ

Table 2

 2 Angles, offsets and amplitudes of the seven lines.

Table 3

 3 Errors on line parameters recovered by the proposed method.

Appendix A. Characterization of the sampling process. Formally, we can write x * φ = (x * φ 1 ) * φ 2 with φ 1 (t 1 , t 2 ) = ϕ 1 (t 1 )δ(t 2 ) and φ 2 (t 1 , t 2 ) = δ(t 1 )ϕ 2 (t 2 ), where ϕ 1 and ϕ 2 are L 1 functions. So, after the first horizontal convolution, using the fact that δ(at) = δ(t)/|a| for any a = 0, we obtain the function:

We can show that, after the second vertical convolution, we get the function b obtained in [START_REF] Blu | Sparse sampling of signal innovations[END_REF]. Figure 15 explains our notations in more details and illustrates the relation between all continuous and discrete variables. In the following we consider a weaker version of Assumption (ii) about ϕ 2 :

Denoting sinc(t 2 ) = sin(πt 2 )/(πt 2 ) the sinus cardinal with sinc(0) = 1, we make the assumption that the discrete filter [START_REF] Starck | Sparse image and signal processing: wavelets, curvelets, morphological diversity[END_REF] h

has compact support of length 2S + 1, for some S ∈ N; that is,

Remark 9. Note that if ϕ 2 is bandlimited, we simply have h[n] = ϕ 2 (n), which reverts back to Assumption (ii).

Let us deduce from these assumptions some other properties satisfied by ϕ 1 and ϕ 2 and their associated discrete filters g and h. First, we have the following proposition:

Proposition 11 (Nyquist-Whittaker-Shannon). The function ϕ 1 , which is periodic and bandlimited, is determined by W degrees of freedom only. That is, with the Altogether, we completely and exactly characterized the sampling process, which involves a continuous blur φ, using the two discrete and finite filters (g[n]) W -1 n=0 and (h[n]) S n=-S . Appendix B. Proof of Proposition 1.

First, we consider the image u obtained by applying the 1-D DFT on every column of u (see [START_REF] Toeplitz | Zur Theorie der quadratischen und bilinearen Formen von unendlichvielen Veränderlichen[END_REF]):

which are the exact Fourier coefficients of the function t → (x * φ 1 )(t, n 2 ), following Assumption (i). Hence, from [START_REF] So | Linear prediction approach for efficient frequency estimation of multiple real sinusoids: algorithms and analyses[END_REF] and

W t dt, we obtain:

Now we apply a 1-D DFT on the first component of the discrete image b (see [START_REF] Tufts | Estimation of frequencies of multiple sinusoids: Making linear prediction perform like maximum likelihood[END_REF]), leading to the elements

Following ( 86) and ( 87), as illustrated in Figure 15, we have

which concludes the proof.

Appendix C. Proof of Proposition 3.

Let us introduce (88) SDP(z) = inf

q∈C N ,q0 0 q 0 : T N (z, q) = T N (q) z z * q 0 0 .

We want to prove that SDP(z) = z A and that the minimum in ( 88) is achieved.

whose components are:

H S ,H S , ∀k = 0, . . . , H S -1 .

Appendix F. Proof of Proposition 9.

First, let us determine the operator norm T M +1 2 = sup

.

By definition, we have z 2

Moreover, we get:

with equality when z = (1, 0, . . . , 0), hence

Let us now decompose the operator T H S as follows:

We directly have T 1 (z)

with equality when q = (1, 0, . . . , 0), hence T 2 2 = H S + 1.

Now, we have

and is achieved when z = 0 and q = (1, 0, . . . , 0), hence

We are now able to derive the operator norms [START_REF] Lassen | Automated Determanation of Crystal Orientations from Electron Backscattering Patterns[END_REF]:

x[m, :]

H , with equality when x = 0, q[0, :] = 0 and q[m, :] = (1, 0, . . . , 0) for all m ∈ 1, M .

Similarly, we have

2 H , with equality when q = 0 and x[:, n 2 ] = (1, 0, . . . , 0) for all n 2 ∈ 0, H S -1 .

We conclude that

= T M +1 2 = M + 1 .

Finally,

, L

Appendix G. Proof of Proposition 7. Let F 1 and F 2 be the mappings:

First, we have

Let us compute the differential function of F :

) and S * 2 (v -P , . . . , v 0 , . . . , v P ) = (v 0 , . . . , v P ). Let us determine the Lipschitz constant of the gradient ∇F :

We are looking for a majoration of each term. We deal with the first one C 1 , the second C 2 being obtained in the same way. Using the inequality (a + b) 2 2a 2 + 2b 2 :

X2 .