Valentin Debarnot

Pierre Weiss

Deep-Blur : Blind Identification and Deblurring with Convolutional Neural Networks

Keywords: Blind deblurring, Blind inverse problems, Identification network, Spatially variant blur, Deep learning, Unrolled network

We propose a neural network architecture and a training procedure to estimate blurring operators and deblur images from a single degraded image. The key assumption is that the forward operators can be parameterized by a lowdimensional vector. The considered models include a description of the point spread function with Zernike polynomials in the pupil plane or product-convolution expansions, which allow to tackle space varying operators. Numerical experiments reveal that the proposed method recovers the blur parameters robustly even for large noise levels. This estimate can then used as an input of an unrolled neural network to deblur the image. It is trained with a specific sampling procedure adapted to a family of parameterized operators. This strategy provides a fast and accurate blind inverse problem solver, requiring nearly no human interaction once the operator parameterization has been set up. It adapts to arbirtrary noise levels and considerably improves alternative blind deblurring softwares in the examples considered in this paper, especially in the high noise regimes.

I. INTRODUCTION

Image deblurring and super-resolution consist in recovering a sharp image x from its blurred and sub-sampled version y = P( Āx), where Ā ∈ R M ×N is a discretized linear integral operator describing the acquisition process and P : R M → R M is some perturbation modelling noise, quantization, saturation... It plays an important role in biomedical and astronomical imaging, where physical phenomena such as diffraction and turbulence strongly reduce the achievable resolution. It also received a continual attention in computer vision, where objects that are moving or out-of-focus create significant artifacts. When the operator Ā describing the optical system is available, this problem can be solved with mature variational inverse problem solvers [START_REF] Chambolle | An introduction to continuous optimization for imaging[END_REF] or data-driven approaches [START_REF] Arridge | Solving inverse problems using data-driven models[END_REF].

However, deriving a precise forward model requires specific calibration procedures, well controlled imaging environments and/or highly qualified staff. In addition, model mismatches result in distorted reconstructions. This can lead to dramatic performance loss, especially for super-resolution applications [START_REF] Alex Von Diezmann | Correcting field-dependent aberrations with nanoscale accuracy in three-dimensional single-molecule localization microscopy[END_REF].

An alternative to a careful calibration step consists in solving the problem blindly: the forward model Ā is estimated together with the sharp image x. Unfortunately, this blind inverse problem is highly degenerate. There is no hope to recover the sharp image without prior assumptions on x and Ā. For instance, assume that Ā is a discrete convolution 1 A preliminary version of this work was published in IEEE ISBI 2021 [START_REF] Debarnot | Deepblur: Blind identification of space variant psf[END_REF].

operator with some kernel h, i.e., y = h x. Then the couple ( h, x) can be recovered only up to a large group of transformations [START_REF] Soulez | Superresolution with opticallymotivated blind deconvolution[END_REF]. For instance, the identity and blurred image are a trivial solution, the image and kernels can be shifted in opposite directions or scaled with inverse factors. Therefore, it is critical to introduce regularization terms both for the operator Ā and the signal x.

The main objective of this work is to design a blind inverse problem solver under the two assumptions below:

• The operator Ā can be parameterized by a lowdimensional vector. In what follows, we let A : R K → R M ×N denote the operator mapping and we assume that Ā = A(γ) for some γ ∈ R K .

• The signal x lives in a family X with some known distribution L X . We propose a specific convolutional neural architecture and a training procedure to recover the couple (γ, x) from the degraded data y and the mapping A(γ). A first network identifies the parameterization γ, while the second uses this parameterization to estimate the image ȳ. This results in an efficient algorithm to sequentially estimate the blur operator and the sharp image x. The network architecture is shown on Fig. 1. Even though the work can be adapted to arbitrary inverse problems beyond image deblurring, we showcase its efficiency for challenging deblurring tasks involving convolutions, but also more advanced space-varying operators.

A. Related works

Solving blind deblurring problems is a challenging task that started being studied in the 1970's [START_REF] Thomas G Stockham | Blind deconvolution through digital signal processing[END_REF]. Fifty years later, it seems impossible to perform an exhaustive review of existing methods and the following description will be lacunary. The prevailing approach is to estimate the original signal and the blur operator by solving variational problems of the form:

inf A∈R M ×N , x∈R N 1 2 Ax -y 2 2 + R A (A) + R x (x), (1)
where R A : R M ×N → R∪{+∞} and R x : R N → R∪{+∞} are regularization terms for the operator and the signal respectively. This problem can then be attacked with various types of alternating minimization procedures [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. The solutions can be understood as maximum a posteriori estimates [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF].

Before the advent of data-driven approaches, the regularizers were carefully designed to target specific features. The point spread functions can be considered as sparse and compactly supported for motion deblurring [START_REF] Krahmer | Blind image deconvolution: Motion blur estimation[END_REF], [START_REF] Fergus | Removing camera shake from a single photograph[END_REF], [START_REF] Dai | Motion from blur[END_REF], [START_REF] Sun | Edge-based blur kernel estimation using patch priors[END_REF], [START_REF] Perrone | Total variation blind deconvolution: The devil is in the details[END_REF], [START_REF] Peyrin | Semi-blind joint superresolution/segmentation of 3d trabecular bone images by a tv box approach[END_REF], [START_REF] Chakrabarti | Analyzing spatially-varying blur[END_REF], [START_REF] Couzinie-Devy | Learning to estimate and remove non-uniform image blur[END_REF]. They are smooth for diffraction limited systems [START_REF] Tony | Total variation blind deconvolution[END_REF] and can also be parameterized with Zernike polynomials in the pupil plane [START_REF] Joseph W Goodman | Introduction to Fourier optics[END_REF], [START_REF] Keuper | Blind deconvolution of widefield fluorescence microscopic data by regularization of the optical transfer function (otf)[END_REF], [START_REF] Aristov | Zola-3d allows flexible 3d localization microscopy over an adjustable axial range[END_REF], [START_REF] Soulez | Blind deconvolution of 3d data in wide field fluorescence microscopy[END_REF], [START_REF] Soulez | Superresolution with opticallymotivated blind deconvolution[END_REF]. The images can be considered as sparse in microscopy and astronomical imaging [START_REF] Mourya | A blind deblurring and image decomposition approach for astronomical image restoration[END_REF], [START_REF] Debarnot | Blind inverse problems with isolated spikes[END_REF] or piecewise constant for natural images. The typical regularizer R x is then the total variation, or more advanced priors on the image gradient [START_REF] Tony | Total variation blind deconvolution[END_REF], [START_REF] Bar | Semi-blind image restoration via mumford-shah regularization[END_REF], [START_REF] Pankajakshan | Blind deconvolution for thin-layered confocal imaging[END_REF], [START_REF] Perrone | Total variation blind deconvolution: The devil is in the details[END_REF], [START_REF] Peyrin | Semi-blind joint superresolution/segmentation of 3d trabecular bone images by a tv box approach[END_REF], [START_REF] Chakrabarti | Analyzing spatially-varying blur[END_REF]. Some authors also advocate for the use of priors on the image spectrul [START_REF] Zachevsky | Blind deblurring of natural stochastic textures using an anisotropic fractal model and phase retrieval algorithm[END_REF], [START_REF] Goldstein | Blur-kernel estimation from spectral irregularities[END_REF], which allow transforming the blind deconvolution problem as a phase retrieval problem under ideal conditions. The most recent variants of these approaches provide excellent results for specific applications such as text deblurring, see e.g., [START_REF] Pan | Deblurring text images via l0-regularized intensity and gradient prior[END_REF]. However, they strongly rely on the detection of specific features (points, edges, textures) which may be absent or inaccurate models of the typical image features. In addition, the problem (1) or its derivatives are usually highly non-convex and the initialization must be chosen very carefully to ensure local convergence to the right minimizer. As a result, these methods need to be carefully tuned for a specific application fields for success.

In the last five years, machine learning approaches have emerged and now seem to outperform carefully handcrafted ones, at least in well controlled conditions. These approaches can be divided in two categories. The first category concerns end-to-end methods that directly estimate the reconstructed image from the observation [START_REF] Christian | Learning to deblur[END_REF], [START_REF] Noroozi | Motion deblurring in the wild[END_REF], [START_REF] Nah | Deep multiscale convolutional neural network for dynamic scene deblurring[END_REF], [START_REF] Aljadaany | Douglas-rachford networks: Learning both the image prior and data fidelity terms for blind image deconvolution[END_REF], [START_REF] Cho | Rethinking coarse-to-fine approach in single image deblurring[END_REF], [START_REF] Mehri | Mprnet: Multi-path residual network for lightweight image super resolution[END_REF]. The second category contains approaches that produce an estimation of the blur operator and the original image. These approaches are specifically tuned for applications in computer vision [START_REF] Sun | Learning a convolutional neural network for non-uniform motion blur removal[END_REF], [START_REF] Christian | Learning to deblur[END_REF], [START_REF] Gong | From motion blur to motion flow: a deep learning solution for removing heterogeneous motion blur[END_REF], [START_REF] Chakrabarti | A neural approach to blind motion deblurring[END_REF], [START_REF] Li | Deep algorithm unrolling for blind image deblurring[END_REF] (motion and out-of-focus blurs) or diffraction limited systems [START_REF] Möckl | Accurate phase retrieval of complex 3d point spread functions with deep residual neural networks[END_REF], [START_REF] Shajkofci | Spatially-variant cnn-based point spread function estimation for blind deconvolution and depth estimation in optical microscopy[END_REF], [START_REF] Shajkofci | Deepfocus: a few-shot microscope slide auto-focus using a sample invariant cnn-based sharpness function[END_REF]. Our work rather falls in the second category since space varying motion blurs are hardly captured by a low-dimensional parameterization. Notice that while the diversity of motion blurs is higher, their inversion is typically easier since their Fourier decay is much slower.

In this list of references, the idea closest to ours is from Shajkofci and Liebling [START_REF] Shajkofci | Spatially-variant cnn-based point spread function estimation for blind deconvolution and depth estimation in optical microscopy[END_REF], [START_REF] Shajkofci | Deepfocus: a few-shot microscope slide auto-focus using a sample invariant cnn-based sharpness function[END_REF]. Therein, the authors estimate a decomposition of the point spread function from a single image using a low dimensional parameterization such as a decomposition over Zernike polynomials. The spatial variations are then estimated by splitting the observation domain in patches where the blur is assumed locally invariant. The image can then be deblurred using a Richardson-Lucy algorithm based on the estimated operator.

B. Contributions

In this work, we propose to use a pair of convolutional neural networks to first estimate the operator parameterization γ ∈ R K and then use this parameterization to estimate the sharp image x ∈ R N with a second convolutional neural network. The first network is the popular ResNet [START_REF] He | Deep residual learning for image recognition[END_REF] as in [START_REF] Shajkofci | Spatially-variant cnn-based point spread function estimation for blind deconvolution and depth estimation in optical microscopy[END_REF]. The second network has the structure of an unrolled algorithm, which offers the advantage of adapting to the forward operator [START_REF] Adler | Solving ill-posed inverse problems using iterative deep neural networks[END_REF], [START_REF] Adler | Learned primal-dual reconstruction[END_REF], [START_REF] Vishal Monga | Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing[END_REF]. We call the resulting algorithm Deep-Blur , see Fig. 1. This work contains various original features:

• It includes space varying blur operators. In particular, we show that this approach allows to use the characterization of an optical system as a low-dimensional subspace of operators proposed in [START_REF] Debarnot | A scalable estimator of sets of integral operators[END_REF], [START_REF] Debarnot | Learning low-dimensional models of microscopes[END_REF]. Most approaches in the literature decompose the observation space into patches and treat each patch independently. In this work, we consider operators with an impulse response that varies continuously in the field of view. • The resulting deblurring network is able to adapt to different forward models and to handle model mismatches naturally. This issue is an important concern for the use of model-based inverse problem solvers [START_REF] Vegard Antun | On instabilities of deep learning in image reconstruction and the potential costs of ai[END_REF], [START_REF] Genzel | Solving inverse problems with deep neural networks-robustness included[END_REF], [START_REF] Ongie | Deep learning techniques for inverse problems in imaging[END_REF].

As will be discussed later, our approach can be seen as an intermediate step between the plug&play algorithms [START_REF] Singanallur V Venkatakrishnan | Plug-and-play priors for model based reconstruction[END_REF], [START_REF] Zhang | Plug-and-play image restoration with deep denoiser prior[END_REF] and the unrolled algorithms [START_REF] Adler | Solving ill-posed inverse problems using iterative deep neural networks[END_REF]. • To train Deep-Blur , we propose a specific sampling procedure of parameterized operators, assuming that the operators parameters live in a compact set. • We evaluate the efficiency, robustness and stability of the proposed approach on various challenging problems, showing that the method is reliable and accurate. The Python implementation of our method is freely available: www.github.com/DeepBlur.

C. Outline

In Section III, we introduce and motivate the mathematical framework. In Section IV, we present the Deep-Blur architecture and learning framework. In Section V, we validate the method on synthetic examples.

II. NOTATION

The bold font refers to vectors, matrices or vector-valued functions. The i-th value of a vector x is denoted x[i] or x i , I M denotes the identity matrix of size M × M . For two operators A and B in R M ×N , we let A, B def.

= Trace(A T B) and the Frobenius norm is defined by A

F = A, A . A family of K operators (A k ) ∈ R M ×N is said to be orthogonal if A k , A k =    0 if k = k 1 otherwise . The Hausdorff distance between two sets X , Y ⊂ R M is defined by dist(X , Y) def. = max sup x∈X inf y∈Y x -y 2 , sup y∈Y inf x∈X x -y 2 .
We let S P -1 = {λ ∈ R P , λ p ≥ 0, P p=1 λ p = 1} denote the unit simplex. Let (γ p ) 1≤p≤P ∈ R K denote a set of P points in R K , the convex cone generated by these vectors is defined by

cone (γ p , 1 ≤ p ≤ P ) def. = P p=1 λ[p]γ p , λ ∈ R P + .
The convex hull is given by

conv (γ p , 1 ≤ p ≤ P ) def. = P p=1 λ[p]γ p , λ ∈ S P -1 .

III. FORWARD MODELLING

In this paper, we assume that the degraded signal y ∈ R M is generated according to the following equation:

y = P(A(γ)x), (2)
where A(γ) : R N → R M is a linear operator depending on the unknown parameter γ ∈ R P . The mapping P : R N → R N can model various deterministic or stochastic perturbations occuring in real systems such as:

• Additive white Gaussian noise: P(ȳ) = ȳ + b with b ∼ N (0, σ 2 I M ). • Poisson noise: In that case, y ∈ N M denotes a number of photons. • Quantization: P(ȳ) = 1 ∆ y • ∆ , where ∆ > 0 is the quantization step.

• Clipping: P(ȳ) = min(ȳ, U ), where U > 0 denotes a saturation parameter indicating the maximum available dynamics of the image. • Combinations: the above perturbations can be combined to produde more accurate models. In this paper, we will focus on the additive white Gaussian noise only, though all these models can be incorporated easily in the proposed framework at the learning stage. A critical aspect of this paper is the parameterization of the forward operator A. We discuss this aspect below.

A. Linear parameterization

The simplest choice is probably to assume that A belongs to a subspace of operators. This leads to the linear model below.

Assumption III.1 (Subspace of operators). The forward operator in Problem (2) belongs to a known K-dimensional subspace of operators:

A(γ) def. = K k=1 γ[k]A k , (3)
with known elementary operators

A k ∈ R M ×N .
Up to a Gram-Schmidt orthogonalization, we can always assume that the operators (A k ) have unit norm and are pairwise orthogonal with respect to the Frobenius scalar product.

The main technical asset of the above model is the following property.

Proposition III.1. Under Assumption III.1, we have

A(γ) F = γ 2 and A(γ) -A(γ ) F = γ -γ 2 for all γ, γ ∈ R P Proof. A(γ) 2 F = K k=1 γ[k]A K k=1 γ[k]A k = K k=1 K k =1 γ[k]γ[k ] A k , A k = K k=1 γ[k] 2 ,
by orthogonality of the family (A k ). The second identity results from the linearity of the parameterization.

To further improve the identifiability of the problem, we will sometimes assume the following:

Assumption III.2. The parameter γ belongs to a known convex cone Γ ⊂ R K , defined with P extreme rays:

Γ = cone(γ p , 1 ≤ p ≤ P ). (4)
Let us now discuss how these assumptions can be used in practice for real imaging systems.

1) Convolution models and eigen-PSF bases: By far the most wide-spread blurring model in imaging is based on convolution operators: the point spread function is identical whatever the position in space. This model is accurate for small fields of view, which are widespread in applications. Assuming that there is no sub-sampling in the model, we can set M = N and Ax = h x for some unknown convolution kernel h.

The convolution model strongly simplifies the blur identification problem since we are now looking for a vector of size N instead of a huge N × N matrix. Yet, the blind deconvolution problem is known to suffer from many degeneracies and possesses a huge number of possible solutions, see e.g. [START_REF] Soulez | Superresolution with opticallymotivated blind deconvolution[END_REF]. To further restrict the space of admissible operators and therefore improve the identifiability, we can expand the kernel h in an eigen-PSF basis. This leads to the following low dimensional model.

Model III.1 (Convolution and eigen-PSFs ). We assume that

A(γ) = K k=1 γ[k]e k x,
where (e k ) is an orthogonal family of convolution kernels coined eigen-PSF basis. This amounts to setting

A k x = e k x
in the frame of Assumption III.1.

In the above model the elementary operators (A k ) are pairwise orthogonal since the basis (e k ) is orthogonal. Defining an eigen-PSF basis can be achieved by computing a principal component analysis of a family of observed or theoretical point spread functions [START_REF] Frisken | Diffraction by a circular aperture as a model for three-dimensional optical microscopy[END_REF]. The authors recently developed the Fiji plugin Eigen-PSF extractor to this end. An example of experimental eigen-PSF basis is shown on Figure 2, top.

Let (h p ) 1≤p≤P define a set of P point spread functions in R N . They can be decomposed over the basis (e k ), and represented by vectors γ p of size K. The set Γ of admissible parameterization can then be defined as the conical hull generated by the projection coefficients (γ p ) as in [START_REF] Alex Von Diezmann | Correcting field-dependent aberrations with nanoscale accuracy in three-dimensional single-molecule localization microscopy[END_REF]. This method nicely captures some physical properties of the point spread functions. For instance, non-negativity and smoothness are preserved by taking convex combinations. Hence, computing a conical hull is a natural way to learn models of impulse responses and interpolate between them.

2) Space variant models and product-convolution expansions: The convolution model III.1 can only capture space invariant impulse responses. When dealing with large field of views, this model becomes inaccurate. One way to overcome this limitation is to use product-convolution expansions [START_REF] Denis | Fast approximations of shift-variant blur[END_REF], [START_REF] Escande | Approximation of integral operators using product-convolution expansions[END_REF], [START_REF] Debarnot | Learning low-dimensional models of microscopes[END_REF], which efficiently encode space-varying systems.

Model III.2 (Product-convolution expansions ).

Let (e i ) 1≤i≤I and (f j ) 1≤j≤J define two orthogonal families of R N . For k = (i, j), we define the k-th elementary operator A k by

A k x = e i * (f j x), ∀x ∈ R N , (5)
where indicates the coordinate-wise (Hadamard) product.

In the above model, the basis (e i ) can still be interpreted as an eigen-PSF basis. Indeed, we have for all locations z ∈ {1, . . . , N }:

A(γ)δ z = I i=1   J j=1 γ[i, j]f j [z]   e i [ • -z].
Hence, we see that each impulse response is expressed in the basis (e i ). The basis (f j ) on its side, can be interpreted as an eigen-space variation basis: it describes how the point spread functions can vary in space. In optical devices such as microscopes, the estimation of the families (e i ) and (f j ) can be accomplished by observing several images of micro-beads [START_REF] Bigot | Estimation of linear operators from scattered impulse responses[END_REF], [START_REF] Debarnot | Learning low-dimensional models of microscopes[END_REF].

An example of experimental product-convolution family is shown in Figure 2 for a wide-field microscope. In that case, the dimension K of the subspace is K = I • J = 16. Again, taking the conical hull of observed decompositions is an efficient way to restrict the family.

B. Nonlinear parameterization and Zernike polynomials

In the previous paragraph, we considered linear parameterizations which have many assets. They are easy to implement and can lead to convex models if one variable (γ or x) has to be estimated while the other is fixed. In addition, they are quite natural models that can be efficiently estimated from experimental data.

However, noise on the data may lead to inaccuracies making it impossible to capture some fine features of the point spread functions. An alternative is to use theoretical ones, leading to the nonlinear world. A popular and effective model in microscopy and astronomy consists in using the Fresnel approximation and to approximate the pupil function with a finite number of Zernike polynomials [START_REF] Hanser | Phase-retrieved pupil functions in wide-field fluorescence microscopy[END_REF], [START_REF] Joseph W Goodman | Introduction to Fourier optics[END_REF]. This model leads to some of the state-of-the-art algorithms for blind deconvolution and super-resolution in microscopy and astronomy [START_REF] Soulez | Blind deconvolution of 3d data in wide field fluorescence microscopy[END_REF], [START_REF] Aristov | Zola-3d allows flexible 3d localization microscopy over an adjustable axial range[END_REF], [START_REF] Sage | Deconvolutionlab2: An open-source software for deconvolution microscopy[END_REF], [START_REF] Nehme | Deepstorm3d: dense 3d localization microscopy and psf design by deep learning[END_REF].

Model III.3 (Fresnel approximation and a Zernike basis).

We assume that the forward model is a convolution with a continuous 3D kernel h(x, y, z). The 3D kernel can be expressed through the 2D pupil function φ as:

h(x, y, z) = φ(w 1 , w 2 ) exp(2iπzd(w 1 , w 2 )) exp(2iπ(w 1 x + w 2 y))dw 1 dw 2 2 with d(w 1 , w 2 ) = (n/λ) 2 -(w 1 + w 2 ) 2 ,
where n is the refractive index of the immersion medium and λ is the wavelength of the observation light.

The complex pupil function φ can be expanded with Zernike polynomials Z k :

φ = exp 2i K k=1 γ[k]Z k ,

PSF basis 97 pixels

Space variations basis 2304 pixels Fig. 2: Examples of eigen-PSF and eigen-space variation bases for a wide-field microscope [START_REF] Debarnot | Learning low-dimensional models of microscopes[END_REF].

where the coefficients γ[k] ∈ R are real numbers2 .

A few examples of slices of point spread functions generated with Model III.3 are displayed in Figure 4. In that case, we used K = 7 Zernike polynomials and set the coefficients γ as uniform random variables with amplitude less than 0.15. As can be seen, a rich variety of impulse responses can be generated with this low dimensional model.

C. The limits of a low-dimensional parameterization

The above models are sufficiently rich to describe most optical devices accurately. In microscopy, they can capture defocus, refractive index mismatches, changes of temperature, tilts of optical components, usual optical abberations with a parameter dimension K smaller than 20.

Notice however that some phenomena can hardly be modelled by low-dimensional parameterizations. In microscopy for instance, diffraction by the sample itself can lead to extremely complicated and diverse forward models better described by nonlinear equations, see e.g., [START_REF] Pham | Three-dimensional optical diffraction tomography with lippmann-schwinger model[END_REF]. Similarly, in computer vision, motion and defocus blurs can vary abruptly with the movement and depth of the objects. The resulting operators would likely require a large number of parameters, which is out of the scope of this paper.

Overall, we see that the proposed contribution is well adapted to the fine characterization of optical systems, but probably do not extend easily to variations induced by the observed objects themselves.

IV. THE PROPOSED METHODOLOGY

In this section, we discuss the core of our contribution: we first detail the neural network architecture and then explain how to train it.

A. The Deep-Blur architecture

We propose to train two different neural networks sequentially: [START_REF] Krahmer | Blind image deconvolution: Motion blur estimation[END_REF]) or textures (power spectrum) followed by adapted inversion procedures. This whole process can be modelled by a set of linear operations (filtering) and nonlinear operations (e.g., thresholding). A convolutional neural network, composed of similar operations, should therefore be expressive enough to estimate the blur parameters. This is the case for the Deep-Blur identification architecture, a ResNet encoder [START_REF] He | Deep residual learning for image recognition[END_REF] shown on Fig. 1, left. It consists of a succession of convolutions, ReLu activation, batch normalization and average pooling layers, allowing to sequentially reduce the image dimensions. The last layer is an adaptive average pooling layer, mapping the output of the penultimate layer to a vector of size K. The weights θ to be trained are the convolution kernels in the convolution layers and the biases. The encoder structure has been proven to be particularly effective for a large panel of signal processing tasks. It is a basic block of the U-Net architecture [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF]. In our experiments, the vector θ has size |θ| = 11, 178, 448.

2) The deblurring network: The proposed deblurring network mimics a Douglas-Rachford algorithm [START_REF] Patrick | Proximal splitting methods in signal processing[END_REF]. It is sometimes called an unrolled or unfolded network. This type of network currently achieves near state-of-the-art performance for a wide range of inverse problems (see e.g., [START_REF] Vishal Monga | Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing[END_REF]). It has the advantages of having a natural interpretation as an approximate solution of a variational problem and naturally adapts to changes of the observation operators.

a) Deep unrolling: For λ > 0, let R γ,λ denote the following regularized inverse:

R γ,λ = (A(γ) T A(γ) + λI) -1 A(γ) T .
For a parameter γ describing the forward operator and an input image y, the Douglas-Rachford algorithm can be described by the following sequence of operations, ran from t = 0 to t = T -1 with T ∈ N.

Algorithm 1 The Douglas-Rachford deblurring network DN

Require: iteration number T ∈ N, operator γ, scale λ ∈ R + z 0 = R γ,λ (y) for all t = 0 → T -1 do x t+1 = PN t (z t ) z t+1 = z t + R γ,λ (2x t -z t ) -x t end for
The initial guess z 0 corresponds to the solution of

z 0 = argmin z∈R N 1 2 A(γ)z -y 2 2 + λ 2 z 2 2 .
It can be evaluated approximately with a conjugate gradient algorithm run for a few iterations (10 in our implementation).

The mapping PN t (ξ t ) : R N → R N can be interpreted as a "proximal neural network". Proximal operators [START_REF] Patrick | Proximal splitting methods in signal processing[END_REF] have been used massively in the last twenty years to regularize inverse problems. A popular example is the soft-thresholding operator, which is known to promote sparse solutions. Here, we propose to learn the regularizer as a neural network denoted PN t , which may change from one iteration to the next. It corresponds to the green layers in Fig. 1.

The parameters ξ that are learned, are the weights ξ t defining the t-th proximal neural network PN t . In our experiments, the networks PN t have the same architecture for all 1 ≤ t ≤ T . We used the current state-of-the-art network used in plug-and-play algorithms called DRUNET [START_REF] Zhang | Plug-and-play image restoration with deep denoiser prior[END_REF], [START_REF] Hurault | Gradient Step Denoiser for convergent Plug-and-Play[END_REF]. We set T = 4 iterations. Each of the 4 proximal networks contain 8, 159, 808 parameters, resulting in a total of |ξ| = 32, 639, 232 parameters to be trained.

b) Differences with plug&play and deep unrolling: The proposed architecture follows closely the usual unrolled algorithms [START_REF] Adler | Solving ill-posed inverse problems using iterative deep neural networks[END_REF], [START_REF] Adler | Learned primal-dual reconstruction[END_REF], [START_REF] Vishal Monga | Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing[END_REF], [START_REF] Li | An algorithm unrolling approach to deep image deblurring[END_REF], [START_REF] Li | Deep algorithm unrolling for blind image deblurring[END_REF]. There is however a major difference: traditionally, these unrolled architectures are trained to invert a single operator. In this paper, we train the network with a family of operators. We recently showed in [START_REF] Gossard | Training adaptive reconstruction networks for inverse problems[END_REF] that this approach only results in marginal performance loss for a given operator if the family is sufficiently small, while providing a massive improvement in adaptivity to all the operators. In a sense, the proposed approach can be seen as an intermediate step between the plug&play priors which are designed to adapt to all possible operators, and the traditional unrolled algorithm adapted to a single one.

B. Designing the objective functions

Now that the neural network architectures have been described, let us turn to the description of the training procedure for the weights θ and ξ.

We propose to first train the identification network IN(θ) alone and then train the deblurring network DN(ξ) with the output of the identification network as an input parameter. This sequential approach presents two advantages:

• First, the memory consumption is lower. The automatic differentiation only needs to store the parameters of the individual networks, instead of both. This significantly reduces the memory footprint. • Second, the identification network can be used independently of the other and it is therefore tempting to train it separately. In metrology applications for instance, where the aim is to follow the state of an optical system through time, the identification network IN is the most relevant brick. In some applications, such as super-resolution from single molecules, the deblurring network could be replaced by a more standard total variation based solver [START_REF] Bredies | Inverse problems in spaces of measures[END_REF], once the operator is estimated. In what follows, we let X ⊂ R N denote a dataset of admissible images/signals and L X denote a sampling distribution over X . We let L Γ denote a sampling distribution supported on the set Γ ⊂ R K . In our experiments, the perturbation P in Equation ( 2) is assumed to be additive white Gaussian noise N (0, σ 2 I M ) with a fixed noise amplitude σ 2 .

1) Training the identification network: Motivated by Proposition III.1, the estimation error for the parameter γ can be measured by the 2 -norm in R K . This leads to the following empirical risk minimization:

inf θ∈R |θ| E x∼L X γ∼LΓ b∼N (0,σ 2 I M ) 1 2 IN(θ)(A(γ)x + b) -γ 2 2 . ( 6
)
Problem ( 6) can be attacked using stochastic gradient descents type algorithms. In our experiments, we used the Adam optimizer [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF] with the default parameters.

2) Training the deblurring network: Once the identification network IN is trained, we turn to the deblurring network. Similarly, we propose to minimize the following empirical risk:

inf ξ∈R |ξ| E x∼L X γ∼LΓ b∼N (0,σ 2 I M ) 1 2 DN(ξ)(y, γ) -x 2 2 , (7)
where y = A(γ)x + b is the degraded image and γ = IN(θ)(y) is the estimated parameter. Of importance, notice that we do not plug the true parameter γ in 7, but rather the estimated one γ. This way, the deblurring network DN can learn to correct model mismatches that may occur at the estimation step.

3) Sampling: Applying the stochastic gradient machinery requires to design the distributions L X and L Γ , i.e., to find a way to draw the images x and the parameters γ.

For L X , we simply use a uniform distribution over a large database of images. Another option could be to simulate the images according to realistic processes for specific applications such as single molecule localization.

C. Sampling a compact set of operators

A nonconventional aspect of this work is the need for a sampling distribution over the parameters in Γ. The easiest way to sample operators is to proceed as with images: take them uniformly at random among a discrete set of observations. Unfortunately, this approach is not always recommendable or even feasible. First, the distribution of observations might be uneven, leading to bad recovery performance for rare operators. In addition, estimating an operator from experimental data can be complicated and time consuming.

To overcome these issues, we first propose to interpolate a few observed operators, to construct a richer set. We then sample the resulting set as evenly as possible with a strategy described hereafter. An example of set was justified in paragraph III-A1 with the convex hull of observations. Alternatives might exist, such as using generative models [START_REF] Zhang | Deblurring by realistic blurring[END_REF] to describe the space of admissible parameters. We do not explore this possibility in this paper.

a) MaxiMin and MaxiSumMin: In what follows, we assume that γ lives in a compact set Γ ⊂ R K . Our objective is to cover the set Γ as uniformly as possible with a fixed budget of L ∈ N points. This problem is tightly related to the sphere packing problem and is known to be extremely hard to solve exactly (see e.g., [START_REF] Horton | Sphere packings, lattices and groups[END_REF], or the 7-th problem of Smale's list for the XXIst century). One possible way to express it is the maximin formulation [START_REF] Mark E Johnson | Minimax and maximin distance designs[END_REF], [START_REF] Pronzato | Minimax and maximin space-filling designs: some properties and methods for construction[END_REF]:

sup ∆=(δ1,...,δ L )∈Γ L min l =l δ l -δ l 2 . ( 8
)
The idea is to spread the points as much as possible while keeping them within the constraint set Γ. The cost function is nonconvex and non-differentiable. It typically possesses a huge number of critical points. As far as we know, there is currently no available software to solve this problem -even approximately -for the dimensions of our problem [START_REF] Pronzato | Minimax and maximin space-filling designs: some properties and methods for construction[END_REF]. In what follows, we propose an original optimization approach scaling well to dimensions of the order D > 100 and a number of sampling points L > 10000. Problem (8) could possibly be solved with a projected subgradient approach. Unfortunately, the convergence would be extremely slow: only the two closest points in the set ∆ would move at each time step. To accelerate the convergence, we propose to solve instead:

sup ∆=(δ1,...,δ L )∈Γ L F (∆) def. = L l=1 min l =l δ l -δ l 2 . ( 9
)
The interest of solving (9) instead of ( 8) can be seen when looking at the expression of the subgradient:

∂ l F (∆) = δ l -δ l * (l) δ l -δ l * (l) 2 , ( 10
)
where l * (l) is the index of the point closest to δ l . The sign = can be replaced by to obtain an element of the Clarke subdifferential. Assuming that a projector Π Γ onto the set Γ is available, we can then attack (9) with as subgradient algorithm:

∆ t+1 = Π Γ L (∆ t + τ t ∂F (∆ t )), (11)
where (τ t ) is a sequence of step-sizes and ∂F (∆ t ) is the vector with coordinates given by [START_REF] Fergus | Removing camera shake from a single photograph[END_REF]. Contrarily to (8), we see that at each iteration, all the points will move simultaneously in a direction opposite to their closest neighbour. The usual theories [START_REF] Boris | Introduction to optimization. optimization software[END_REF] for subgradient descents recommend a decaying step-size such as τ t = c √ t . In this work, we experimented various choices (e.g., Barzilai-Borwein) and finally opted for a constant step-size for this nonconvex, non differentiable function. Such a choice does not lead to a convergent sequence (∆ t ) in general. It typically finishes oscillating around a local maximizer of the energy with a radius that depends on the step-size. This is sufficient for our application since we -in any cases -end up close to a local maximizer only and do not require a precise solution.

b) Handling convex hulls: In the specific case where Γ = conv(γ p , 1 ≤ p ≤ P ), we have Γ = ΓS P -1 , where Γ = [γ 1 , . . . , γ P ] is the matrix formed by concatenating all the vectors γ p . Computing the orthogonal projection on Γ is nontrivial. To solve this issue, we propose to replace (9) by sup Λ=(λ1,...,λ L )∈S L P -1

F (ΓΛ), (12)
and use a subgradient descent as well:

∆ k+1 = Π S L P -1 (Λ k + τ k Γ T ∂F (ΓΛ k )). ( 13
)
The projection on the simplex Π S P -1 can be evaluated exactly in O(P ) operations using dynamic programming [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF]. The algorithm is stopped when the minimal distance between all pairs of points do not vary significantly after 10 iterations.

Letting Λ = ( λ 1 , . . . , λ L ) denote the last iterate of the algorithm, the points δ l in the previous formulation are given by δ l = Γ λ l .

c) Evaluating the subgradient: The most complicated part of the code is to find the closest neighbor of each point δ l . A naive implementation would result in a complexity O(L 2 ), which is not viable for large L (say larger than 10, 000). Hopefully, finding nearest neighbors efficiently is a problem of high importance and a very efficient implementation was proposed in [START_REF] Johnson | Billion-scale similarity search with gpus[END_REF]. The authors propose a CPU/GPU implementation returning approximate or exact nearest neighbors scaling up to hundreds of billion vectors in high dimension.

d) Examples: To conclude this section, we provide a simple example of the algorithm dynamics in Fig. 3

V. NUMERICAL EXPERIMENTS

In this section, we illustrate the different ideas proposed in this paper. In all our experiments, we trained the neural networks using the MS Coco dataset [START_REF] Lin | Microsoft coco: Common objects in context[END_REF]. It contains 118,287 images in the training set and 40,670 images in the test set.

A. The interest of operator space sampling

We illustrate the importance of carefully sampling the set of admissible operators using the method described in Section IV-C. To this end, we use the blur operator generated by Model III.2. The functions (e i ) 1≤i≤I and (f j ) 1≤j≤J are generated from a real wide-field fluorescence microscope. The dimension of the subspace is K = 4 × 4 and the cone Γ is defined by 400 operators. Details about this construction can be found in [START_REF] Debarnot | Learning low-dimensional models of microscopes[END_REF]. We use Equation [START_REF] Perrone | Total variation blind deconvolution: The devil is in the details[END_REF] to generate a set of 10,000 well spread admissible operators. The coefficients γ are then drawn uniformly at random among this set. We also add Gaussian noise with amplitude σ = 0.01.

We first train the identification network IN with operators among the discretization of the cone Γ containing 10,000 points. We train a second network with a more naive approach where IN is trained using only the 400 blur operators defining the set Γ. After training over 150,000 iterations with batches of 16 images, we evaluate the results on 100 images and 100 operators randomly sampled from the discretization of Γ, bringing the total number of test samples to 10,000. We report the SNR between the true operators and the estimations by the network IN in Table I. Sampling from the discretization of the set Γ comes with almost no additional computational cost and allows to significantly increase the SNR of the estimated blur operators.

B. Convolution operators

In this section, we evaluate the accuracy of the identification and deblurring networks for convolution (i.e. space invariant) operators. We assess them for images generated with point spread functions expanded in Zernike polynomial.

1) Identifying convolution operators: In this paragraph, we assess the ability of a residual network to identify the point spread function generated by Model III.3. A similar study was carried out in [START_REF] Shajkofci | Deepfocus: a few-shot microscope slide auto-focus using a sample invariant cnn-based sharpness function[END_REF] with K = 3 coefficients. Here, we extend the study to K = 7 coefficients allowing to represent the following aberrations in the Noll nomenclature [START_REF] Robert | Zernike polynomials and atmospheric turbulence[END_REF]: defocus, primary astigmatism, primary coma, trefoil and primary spherical.

We generate random PSFs by drawing the coefficients γ[k] (see Model III.3) uniformly in the range [-η, η]. The higher η, the more spread and oscillating the PSF. Hence η can be interpreted as a measure of PSF complexity. The model was trained for a value of η = 0.15.

Figure 4 shows the identification results for 3 images taken at random from the test set and 3 operators taken at random in the operator set. On these examples, the network provides faithful estimates despite an important noise level and images with little contents.

To further characterize the network efficiency, we measure the distribution of signal-to-noise-ratio (SNR) in the noiseless regime. For a kernel h, the error of the estimated kernel ĥ is defined by: SNR(h, ĥ) = -10 log 10 ĥ -

h 2 2 h 2 2 . (14)
Fig. 5 summarizes the conclusions. In average, the identification network outputs estimates with a relative error below 5%.

Finally, we study the stability to the noise level σ in Fig. 6a -and to the PSF complexity η in Fig. 6b. As can be seen, the identification outputs predictions with less than 10% error with probability larger than 0.5 up to a large noise level of σ = 0.1. The dependency on the kernel's complexity, measured through the Zernike polynomials amplitude η is very clear with typical errors below 2% for η < 0.1 and then a relatively fast increase. It is nonetheless remarkable that the identification returns estimates with less than 15% error for η = 0.2, which produces significantly more complex PSFs than those observed during the training phase, showing some ability of the network to extrapolate.

2) Evaluating the deblurring network: In this section, we evaluate the performance of the proposed deblurring network for convolution operators defined using the Fresnel approximation. Figures 7,8, 9 display some deconvolution results for different methods. The corresponding image quality measures are displayed in Tables II, III, IV.

Notice that this problem is particularly involved: there is complete loss of information in the high frequencies since the convolution kernels are bandlimited and we treat different noise levels up to rather high values (here σ = 0.05 for images in the range [0, 1]). Despite this challenging setting, it can be seen both perceptually and from the SSIM (Structural Similarity Index Measure) that the image quality is significantly improved whatever the noise level. It is also remarkable to observe that the proposed network architecture allows to treat images with different noise levels. This is an important feature of the DRUNET used as a proximal network [START_REF] Zhang | Plug-and-play image restoration with deep denoiser prior[END_REF].

We also propose some comparisons with other methods from the literature. Whenever possible, we optimized the hyperparameters by hand for each noise level to produce the best possible output. We chose the following methods:

• The 0 -gradient prior [START_REF] Pan | Deblurring text images via l0-regularized intensity and gradient prior[END_REF], [START_REF] Pan | Blind image deblurring using dark channel prior[END_REF]. This method is one of the state-of-the-art handcrafted blind deblurring method. An efficient implementation was recently proposed in [START_REF] Anger | Blind Image Deblurring using the l0 Gradient Prior[END_REF].

• In [START_REF] Goldstein | Blur-kernel estimation from spectral irregularities[END_REF], the authors proposed a kernel estimation method based on the assumption that the image spectrum amplitude has a specific decaying distribution in the Fourier plane. The kernel estimation then boils down to a phase retrieval problem. An efficient implementation was recently proposed in [START_REF] Anger | Estimating an Image's Blur Kernel Using Natural Image Statistics, and Deblurring it: An Analysis of the Goldstein-Fattal Method[END_REF].

• In [START_REF] Anger | Efficient blind deblurring under high noise levels[END_REF], the authors propose a regularized version of the 0 -gradient prior allowing to treat blurred images containing higher noise levels. This method also implements the FFDNET [START_REF] Zhang | Ffdnet: Toward a fast and flexible solution for cnn-based image denoising[END_REF] as a pre-processing step to remove the noise prior to the estimation. In this table we evaluated the identification network for 100 images (left to right) and 100 kernels (top to bottom) with no noise σ = 0. As can be seen, there are horizontal and vertical stripes. This means that some images and some kernels make the identification problem easier or harder. In the middle: an image making the identification problem hard (column 23). On the right: a kernel making the identification harder (row 65).

0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 • We also tested two state-of-the art neural network approaches. The first one is the current leader (as of May 15, 2022) of the Go-Pro deblurring challenge called NAFNET [START_REF] Chen | Simple baselines for image restoration[END_REF].

• The second one is the Deblur-Gan V2 [START_REF] Orest Kupyn | Deblurgan: Blind motion deblurring using conditional adversarial networks[END_REF], [START_REF] Orest Kupyn | Deblurgan-v2: Deblurring (orders-of-magnitude) faster and better[END_REF]. As can be seen from the tables and the perceptual results, our method (Deep-Blur ) outperforms the other ones by a really significant margin. This is particularly true for the deep learning approaches, which are unable to improve the sharpness significantly and even fail completely for high noise levels. The comparisons are somewhat unfair since we test our method for the conditions it was trained for (diffraction blurs and grayscale images), while the two other deep learning methods were trained on motion blurs and color images. This probably illustrates the lack of robustness of neural network when using conditions escaping the ones used during the training approach. The handcrafted methods (Goldstein-Fattal and 0 gradient prior) significantly improve the image quality in the low noise regime, but fail to do so for high noise levels. a) Training on the true or estimated operators?: At training time, we can feed the unrolled deblurring network with the operator that was used to synthetize the blurry image, or the one estimated using the identification network. The potential advantage of the second option is to train the proximal networks correcting model mismatches. We tested both solutions on two different operator families. It turns out that they led to near undistinguishable results overall in average. The most likely explanation for this phenomenon is that the model mismatches produced by the identification network cannot be corrected with the proximal networks.

b) The interest of training on a family: The proposed deblurring algorithm can be seen as an intermediate step between the plug&play priors which are designed to solve arbitrary inverse problems and the traditional unrolled networks, which are designed for a specific operator. In this paragraph, we investigate the merits of these different approaches empirically. This study is closely related to the one in the preprint [START_REF] Gossard | Training adaptive reconstruction networks for inverse problems[END_REF] for a different application.

To this end, we construct a plug&play deblurring algorithm using the Douglas-Rachford structure described in Algorithm 1. The proximal networks are all identical and correspond to the state-of-the-art DRUNET denoisers used in plug&play approaches [START_REF] Zhang | Plug-and-play image restoration with deep denoiser prior[END_REF]. It was trained specifically for denoising white Gaussian noise. The parameters of the algorithm are the number of iterations, the scaling parameter λ and a standard deviation for the noise level σ. Those parameters were tuned by hand to obtain the best possible signal-to-noise ratio. This optimum is attained with 8 iterations, λ = 10 -1 and σ = 10 -1 .

We also trained our deblurring network DN for a single convolution operator with a Airy pattern. This pattern is obtained by setting γ = 0 in Model III.3.

We then test each deblurring network for 100 differents pairs of images and operators and display the average gain and standard deviation in Table V.

All 3 deblurring networks depend on the operator parameterization γ. We evaluate the quality of the networks independently of the possible model mismatches and therefore feed all the networks with the true parameter used to synthetize the blurry images.

As can be seen in Table V, the plug&play approach takes significantly more time than the lighter deblurring networks we chose for this paper. The reason is twofold: first we use 4 iterations in our deblurring networks instead of 8 for the plug&play approach. Second, the DRUNET used in [START_REF] Zhang | Gaussian approximations of fluorescence microscope point-spread function models[END_REF] contains significantly more parameters than the one used in our implementation. In both cases, memory issues justify this choice. We could likely increase the image quality of Deep-Blur further using more iterations and more parameters in the DRUNET, but this would come at a higher memory/energy cost.

Second, the deblurring network trained specifically on a Airy pattern outperforms the two other methods when tested on a Airy pattern, which is to be expected. However, the gain is anecdotic compared to the same deblurring network trained on a family of operators. In contrast, the deblurring network trained on a family of operators signicantly outperforms the one trained on the Airy disk when tested on multiple operators. This study therefore confirms the findings in [START_REF] Gossard | Training adaptive reconstruction networks for inverse problems[END_REF], where it was shown that the price of adaptivity to different operators was affordable. The plug&play approach on its side performs remarkably well for a "universal" method. Its main downsides are the computing times as well as a slight performance decrease.

C. Product-convolution operators

To finish the numerical experiments, we illustrate how the proposed ideas perform on product-convolution operators. We first illustrate the performance of the identification network. We trained the identification network on natural images from the MS Coco dataset, but evaluate it on biological images from microscopes. We selected 6 images: Image 1 is an histopathology of angiolipoma [START_REF]File:histopathology of angiolipoma.jpg -wikimedia commons, the free media repository[END_REF], Image 2 is an histopathology of reactive gastropathy [START_REF]File:histopathology of reactive gastropathy.jpgwikimedia commons, the free media repository[END_REF], Image 3&4 are actin filaments within a cell [START_REF]File:std depth coded stack phallodin stained actin filaments.png -wikimedia commons, the free media repository[END_REF], Image 5 is an slice of a spheroid from [START_REF] Lorenzo | Live cell division dynamics monitoring in 3d large spheroid tumor models using light sheet microscopy[END_REF] and Image 6 is a crop of a podosome obtain on a wide-field microscope [START_REF] Bouissou | Podosome force generation machinery: a local balance between protrusion at the core and traction at the ring[END_REF].

The blur operators are generated by Model III.2 using K = 16 parameters. This model was estimated using a real optical system [START_REF] Debarnot | Learning low-dimensional models of microscopes[END_REF]. Fig. 10 shows the identification results. The blur coefficients predicted by the Deep-Blur identification are accurate estimates in all cases. In average, the SNR is much higher than in the previous experiment, which can likely be explained by a smaller dimensionality of the operators' family. In all cases, the image quality is significantly improved despite an additive white Gaussian noise with σ = 1 • 10 -2 . This is remarkable since this type of images is really different from the typical computer vision images found in the MS COCO dataset.

VI. CONCLUSION

We proposed a specific neural network architecture to solve arbitrary blind inverse problems. We evaluated its performance carefully on blind deblurring problems with space invariant and space varying operators. A key assumption is to have access to a forward model that depends on a set of parameters. The network first estimates the unknown parameters describing the forward model from the measurements with a ResNet architecture. In a second step, an unrolled algorithm solves the inverse problem with a forward model that was estimated at the previous step. After designing a careful training procedure, we showed a significant advantage of the proposed approach in terms of robustness to noise levels and adaptivity to a vast family of operators and conditions not seen during the training phase. The deblurring network can be seen as an intermediate step between a plug&play algorithm which solves any inverse problems and the traditional unrolled algorithms which solve a specific one.

Fig. 1 :

 1 Fig. 1: The Deep-Blur architecture. The first part of the network identifies the parameterization γ. In this paper, we use a ResNet architecture. The estimated parameterization γ is given as an input of a second deblurring network. This one is an unrolled Douglas-Rachford algorithm. The yellow blocks are convolution layers with ReLu and batch normalization. The red ones are average pooling layers. The green ones are regularized inverse layers of the form xt+1 = (A * (γ)A(γ) + λI) -1 A(γ)y. The violet blocks are UNET like neural networks with weights learned to provide a sharp image x.

•

  The first one denoted IN is an identification network. It depends on weights θ. The mapping IN(θ) : R M → R K takes as an input a degraded image y ∈ R M and provides an estimate γ of γ in R K . • The second one denoted DN is a deblurring network. It depends on weights ξ. The mapping DN(ξ) : (R M , R K ) → R N takes as input parameters the blurry image y and an operator parameterization γ. It outputs an estimate x of the sharp image x. 1) The identification network: Traditional estimation of a blur kernel relies on the detection of cues in the image such as points (direct observation), edges in different orientations (Radon transform of the kernel

  in 2D.

Fig. 3 :

 3 Fig.3: Dynamics of the maxisummin algorithm for L = 500 points. The pair of numbers below the figures are the iterate number and the distance between the closest pair of points in the design. As can be seen, after 10,000 iterations, the points are near perfectly spread in the convex hull of 8 points. Here, we set τt = 10 -3 .

Fig. 4 :Fig. 5 :

 45 Fig.4: Examples of results for the identification network with convolution kernels defined through Fresnel approximation. Top: the original and blurred and noisy 400x400 images. Bottom: the true 31x31 kernel used to generate the blurry image and the corresponding estimation by the neural network. Notice that there is a significant amount of white Gaussian noise added to the blurred image. The image boundaries have been discarded from the estimation process to prevent the neural network from using information that would not be present in real images.

Fig. 6 :

 6 Fig. 6: Stability of the kernel estimation with respect to noise level (left) and amplitude of the Zernike coefficients in the noiseless regime(right).

Fig. 7 :Fig. 8 :Fig. 9 :

 789 Fig. 7: Deep-Blur in action in the noiseless setting. Quantitative evaluations are reported in TableII. When available, the estimated blur kernel is displayed at the bottom-right. First row: original images. Second row: blurry-noisy images. Third row: Deep-Blur. Fourth row:[START_REF] Anger | Blind Image Deblurring using the l0 Gradient Prior[END_REF] Fiveth row:[START_REF] Anger | Estimating an Image's Blur Kernel Using Natural Image Statistics, and Deblurring it: An Analysis of the Goldstein-Fattal Method[END_REF] Sixth row:[START_REF] Anger | Efficient blind deblurring under high noise levels[END_REF] Seventh row:[START_REF] Chen | Simple baselines for image restoration[END_REF] Eith row:[START_REF] Orest Kupyn | Deblurgan-v2: Deblurring (orders-of-magnitude) faster and better[END_REF]

Fig. 10 :

 10 Fig. 10: Deep-Blur applied to spatially varying blur operators on microscopy images (not seen during training). The blur operators are sampled from a family estimated using a real wide-field microscope, as described in Paragraph V-A. First row: the original images. Second row: blurry-noisy images. Third row: the blind deblurring result with Deep-Blur . The SSIM of the resulting deblurred image is displayed below. Fourth row: The true blur operator. We display 4 evenly spaced impulse responses in the field of view. Fifth row: The estimated blur operator. The SNR of the estimated kernel is displayed in the caption in dB.

TABLE I :

 I Importance of sampling the cone Γ of admissible operators. SNR between the true operator and the one estimated by IN for two different set of operators used during training.

TABLE II :

 II Reconstruction results in the noiseless setting σ = 0.

				im 1	im 2	im 3	im 4	im 5	im 6
			y	0.47	0.43	0.40	0.50	0.46	0.49
			ours	0.58	0.84	0.67	0.81	0.79	0.81
			[77]	0.53	0.63	0.53	0.52	0.53	0.64
	SSIM	x	[78]	0.55	0.54	0.53	0.71	0.55	0.58
			[79]	0.52	0.52	0.43	0.54	0.50	0.48
			[80]	0.50	0.48	0.45	0.54	0.49	0.52
			[81]	0.52	0.42	0.38	0.39	0.43	0.56
			ours	8.66	14.48 13.38 11.62	12.19 15.80
	SNR	ĥ	[77] [78]	-5.28 0.94	-0.25 1.19	-1.08 0.59	-2.61 0.20	1.83 0.83	2.60 2.25
			[79]	-0.28	1.37	0.57	-0.36	3.57	2.85

TABLE III :

 III Reconstruction results σ = 10 -2 .

				im 1	im 2	im 3	im 4	im 5	im 6
			y	0.29	0.25	0.25	0.29	0.31	0.36
			ours	0.52	0.50	0.52	0.72	0.60	0.58
			[77]	0.47	0.48	0.41	0.49	0.53	0.39
	SSIM	x	[78]	0.36	0.36	0.32	0.41	0.40	0.36
			[79]	0.49	0.53	0.38	0.56	0.53	0.47
			[80]	0.0	0.32	0.0	0.0	0.0	0.0
			[81]	0.32	0.27	0.25	0.30	0.31	0.38
			ours	7.05	15.00 12.88 11.03 13.01 19.63
	SNR	ĥ	[77] [78]	-7.66 0.32	-4.85 0.42	-3.01 0.38	-6.81 0.50	-0.21 0.71	-3.50 0.62
			[79]	-1.63	-0.53	0.94	-2.83	2.50	2.83

TABLE IV :

 IV Reconstruction results σ = 5 • 10 -2 .

TABLE V :

 V Average SSIM gain for the different networks. The gain was evaluated on 100 random test images.

In this model, only the phase of the pupil function varies. In all generality, the amplitude could vary as well at a slower rate. Most models in the literature assume a constant amplitude.

ACKNOWLEDGMENTS

The authors wish to thank Emmanuel Soubies and Thomas Mangeat for fruitful discussions. This work was supported by the ANR Micro-Blind and by the ANR LabEx CIMI (grant ANR-11-LABX-0040) within the French State Programme "Investissements d'Avenir". The authors acknowledge the support of AI Interdisciplinary Institute ANITI funding, through the French "Investing for the Future-PIA3" program under the Grant Agreement ANR-19-PI3A-0004. This work was performed using HPC resources from GENCI-IDRIS (Grant 2021-AD011012210R1).