<figure><img src="image_1.png" /></figure>

HAL
open science

Learning low-dimensional models of microscopes

Valentin Debarnot, Paul Escande, Thomas Mangeat, Pierre Weiss

► To cite this version:

Valentin Debarnot, Paul Escande, Thomas Mangeat, Pierre Weiss. Learning low-dimensional models of microscopes. IEEE Transactions on Computational Imaging, inPress, 7, pp.178-190. 10.1109/TCI.2020.3048295. hal-02445642v3

HAL Id: hal-02445642

https://hal.science/hal-02445642v3

Submitted on 15 Feb 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Learning low-dimensional models of microscopes

Valentin Debarnot, Paul Escande, Thomas Mangeat, Pierre Weiss

***Abstract***—We propose accurate and computationally efficient procedures to calibrate fluorescence microscopes from micro-beads images. The designed algorithms present many original features. First, they allow to estimate space-varying blurs, which is a critical feature for large fields of views. Second, we propose a novel approach for calibration: instead of describing an optical system through a single operator, we suggest to vary the imaging conditions (temperature, focus, active elements) to get indirect observations of its different states. Our algorithms then allow to represent the microscope responses as a low-dimensional convex set of operators. This approach is deemed as an essential step towards the effective resolution of blind inverse problems. We illustrate the potential of the methodology by designing a procedure for blind image deblurring of point sources and show a massive improvement compared to alternative deblurring approaches both on synthetic and real data.

***Index Terms***—Space varying PSF, blind deblurring & deconvolution, machine learning, operator modelling, product-convolution, super-resolution, sparsity

## I. INTRODUCTION

Many recent breakthroughs in optics pertain to the field of computational microscopy: computers play a critical role to generate images. This evolution allowed to observe objects with unprecedented contrasts, temporal/spatial resolutions or gave access to new quantitative features. To name a few examples, let us mention Single Molecule Localization Microscopy (SMLM), Structured Illumination Microscopy (SIM), Total Internal Reflection Fluorescence microscopy (TIRF) or Stimulated Emission Depletion (STED) microscopy.

A common prerequisite for these techniques is the design of an accurate mathematical model of the optical system. This step is critical since the generation of images usually relies on an explicit or implicit inversion of this model. The advent of these microscopes therefore makes it more and more important to finely characterize their transfer function.

By far, the dominant models in the image processing literature are space invariant systems: the point spread function (PSF) is identical wherever in space. While simplifying the theoretical, numerical and experimental aspects, this assumption is however often unrealistic. Following [1], the space invariance is approximately valid only under very restrictive assumptions. There is clear theoretical and experimental evidence that the variations of the PSF need to be taken into account along the optical axis [2], [3] and the lateral axis for large numerical apertures (see for instance Fig. 3a). Neglecting this aspect can have dramatic consequences. For

instance, it was shown in [4], [5] that this approximation can severely damage the reconstruction of images in single molecule localization imaging, with localization errors of more than 20% for a displacement of less than 200nm. The effects would be even more stringent for large fields of views which are a current challenge with the improved quality of sCMOS detectors. These model mismatches can significantly downgrade the performance of all computational microscopy systems and it is hence critical to finely estimate the optical response of the system.

a) *Existing works:* A well spread approach to describe the response of an imaging system consists in using Fourier optics [6]. This theory provides a nice description of the system through the pupil function of the objective. In this domain, it is possible to derive mathematical models of space variant PSFs [5], [7], [8] and to infer the parameters of these models (e.g. Zernike coefficients) from experimental data. There are however three limitations to these approaches. First, they are often based on parameters such as the numerical aperture, the wavelengths, or many other physical quantities and the models become significantly more complex for large numerical apertures [1], [9]. The more parameters, the more precise the model, but the harder it becomes to finely characterize them experimentally. In addition, some active components such as micro-mirrors introduce additional perturbations which cannot be easily modeled or inferred. The second problem comes from numerical considerations: the dependency between the model and its parameters is non-linear, which inevitably leads to non-convex estimation procedures, leading to local minima and additional inaccuracies. Finally, it is important to notice that what matters is not just the PSF related to the optical elements, but rather the so-called effective PSF [10], which comprises other effects, such as pixel integration and thermal diffusion across sensors.

Instead of relying on a physical model it is possible to directly estimate the optical response from the data. We will follow this approach in this paper. Imaging fluorescent micro-beads in a cover-slide in 2D or a cylinder of agarose in 3D gives a partial idea of the system by providing an access to a few sampled and noisy scattered PSFs. This information can be used to estimate a space invariant system by averaging multiple micro-beads [11]–[13]. When the images are aliased it is even possible to obtain a super-resolved estimation [14], [15]. It becomes more delicate to estimate a space variant system. A few researchers - especially in the field of astrophysics <fnref n="1" /> - have addressed this issue. The general

V. Debarnot & P. Weiss are with IMT, CNRS and Université de Toulouse, France. P. Escande is with I2M, CNRS and Université de Aix-Marseille, and T. Mangeat is with CBI, CNRS and Université de Toulouse

<footnote n="1">¹In astrophysics, the PSFs variations may be due to weak gravitational lensing and reveal distant massive galaxies.</footnote>
framework is the following: a parameterized PSF model is designed either from physics equations or from the data itself. The observed PSFs are then interpolated to cover the whole field of view. In [16], [17], the authors propose to decompose the PSFs over a low dimensional basis and to interpolate the coefficients using thin-plate splines. A subset of the authors proved that this method was minimax optimal in [18] and we will propose a refined version in this paper. It is also possible to use more advanced interpolation methods, using matrix factorization techniques [19], which share similarities with the proposed approach or optimal transport [15], but this method would not scale computationally to the large field of views considered in here.

b) Our contribution: While there now exists a solid theoretical and algorithmic framework to estimate space varying optical responses of optical systems, these methods are often tested on synthetic data that do not reproduce all the complexity of real microscopy images. For instance, the estimation of a PSF requires a very careful treatment of the background and of the noise. Its interpolation requires specific care to avoid obtaining unrealistic results far from the observed responses. Our first objective is to provide *precise estimation algorithms adapted to real data emanating from fluorescence microscopy*. In particular, we propose an algorithm to promote realistic PSFs encoding properties such as nonnegativity or smoothness properties by estimating a conical hull of projection coefficients. This approach shares similarities with [19], but differs in that the features of realistic PSFs are directly learned from the data rather than defined as priors. The second and arguably most original contribution is a *new way to calibrate an optical system by learning all its possible states*. Instead of trying to estimate a single operator to describe the microscope, we propose to learn a whole family of possible states by varying the experimental conditions, following the recent theoretical work [20]. We finally show how the proposed methodology allows to *design new algorithms for blind image deblurring<fnref n="2" /> of point sources*, which significantly outperform more traditional models proposed in the literature.

## II. OPERATOR ESTIMATION

### A. Notation

In all the paper, bold fonts refer to vectors, matrices or vectorial functions while regular fonts refer to scalar numbers or functions. The *i*-th value of a vector *x* is denoted either $x_i$ or $x[i]$. The $\ell^p$ norm of a vector *x* is denoted $\|x\|_p$. The value of a function *f* is $f(x)$ and its $\ell^p$-norm is denoted $\|f\|_p$. The delta Dirac function at a position $x \in \mathbb{R}^d$ is denoted $\delta_x$. In all the paper, the integers *I*, *J*, *K*, *M* and *N* refer to a number of components described in Table I.

<table><caption>TABLE I: Notation.</caption><thead><tr><th>Symbol</th><th>Meaning</th></tr></thead><tbody><tr><td>I</td><td>size of PSF basis</td></tr><tr><td>J</td><td>size of space variations basis</td></tr><tr><td>K</td><td>number of observed microbeads images</td></tr><tr><td>M</td><td>number of observed microbeads</td></tr><tr><td>N</td><td>number of pixels of an image</td></tr></tbody></table>

### B. Preliminaries

A space varying blurring operator $H: L^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)$ can be seen as a linear integral operator mapping an image $u$ to the degraded image $Hu$ through the formula

$$Hu(\mathbf{x}) = \int_{\mathbb{R}^d} L(\mathbf{x}, \mathbf{y})u(\mathbf{y}) \, d\mathbf{y}. \quad (1)$$

The function $L(\cdot, \cdot)$ is called *kernel* of the operator. It describes the impulse response of the system at every location $z \in \mathbb{R}^d$ of the image domain since:

$$(H\delta_z)(\cdot) = L(\cdot, z). \quad (2)$$

The PSF $S(\cdot, z)$ of the system at $z$ is defined as the impulse response centered at 0, i.e.

$$S(\cdot, z) = L(\cdot - z, z). \quad (3)$$

The function $S$ is called the space varying impulse response of the system. This work is based on two important assumptions.

**Assumption 1 (PSF approximation).** Every PSF in the field of view is well approximated by its projection over a low-dimensional orthogonal basis $(h_i)_{1 \le i \le I}$, i.e.

$$S(\cdot, z) \approx \sum_{i=1}^{I} \langle S(\cdot, z), h_i \rangle h_i. \quad (4)$$

This assumption is valid both from a theoretical and an empirical viewpoint. It is indeed well known that any smooth function can be well approximated by its projection on a low dimensional subspace. Typical bases include splines or low frequency Fourier atoms [21]. In practice, we can also construct the family $(h_i)$ by computing the principal component analysis of a family of sampled PSFs. The numerical experiments performed in the paper reveal that for our imaging systems, as little as $I = 5$ elements are enough to capture all possible PSFs accurately. In addition, restricting the PSFs to live on a low dimensional subspace is an efficient method to denoise them as will be illustrated in the numerical section.

**Assumption 2 (PSF variations).** Letting

$$\alpha_i(\mathbf{x}) := \langle S(\cdot, \mathbf{x}), h_i \rangle \quad (5)$$

denote the *i*-th coefficient of the PSF at $\mathbf{x} \in \mathbb{R}^d$, we assume that $\alpha_i$ varies slowly in space.

This hypothesis means that the PSFs vary smoothly in space. It can be substantiated experimentally using arrays of micro-beads for instance, see e.g. [4], [5], [22].

<footnote n="2">²Throughout the paper, we will prefer the term *deblurring* to the more common *deconvolution*. We believe that the first one should be preferred whenever referring to space varying systems: by definition a convolution is space invariant and would be improper in our context.</footnote>
Under Assumptions 1 and 2, a space varying operator is completely characterized by the pair $(\alpha_i, h_i)_{1\le i\le I}$. Of interest, this representation of the operator also leads to fast numerical computations using a structure called product-convolution.

**Proposition 1** (Product-convolution [23]–[26]). Assume that a blurring operator $H$ has a space varying impulse response $S$ defined by the tensor:

$$S(\mathbf{x}, \mathbf{y}) = \sum_{i=1}^{I} \alpha_i(\mathbf{x})h_i(\mathbf{y}),$$

then for any image $u$, we have

$$Hu = \sum_{i=1}^{I} h_i \star (\alpha_i u),$$

where the symbol $\star$ stands for the convolution operator.

Hence, the numerical complexity of computing a space-varying operator is just $I$ times the one of a convolution, which can be achieved efficiently using Fast Fourier Transforms for instance in a way independent of the PSF size.

To the best of our knowledge, this decomposition was first introduced in [23] and has been revisited and improved for the last two decades. It can be found under other names such as deformable filtering, filter flow, parallel product-convolution, anisoplanatic convolutions... Recent works [25], [26] provide an overview of existing choices of generic elementary functions $(\alpha_i)$ and $(h_i)$. The precise approximation rates of product-convolution operators have been studied in [25] and we refer to the previous references for more insight on these structures. The product-convolution expansion has been used to restore images when the blur operator is known [2], [27]–[29], but less attention has been paid to the estimation of these operators. In [30], the authors assume that the collection of elementary PSFs $(h_i)$ is known. Then, they estimate the space variations and the original signal. The approach of this paper is rather different since we aim at estimating both families of space variations and PSFs beforehand. In this work, we follow the theory developed in [20] and learn these functions from the acquisitions.

### C. Estimating a single operator

Under Assumptions 1 and 2, the problem of estimating the operator reduces to recovering the low dimensional bases $(h_i)$ and $(\alpha_i)$, or at least their discretized counterparts $(h_i)$ and $(\alpha_i)$. In this paragraph, we describe the general principle of the estimation of a blurring operator from a single image of fluorescent micro-beads. The general process is described in Fig. 1.

The first step consists in extracting the most relevant PSFs in the form of small patches (see Appendix A). Then the background is estimated and removed on each patch independently to avoid biases in the PSFs estimation. We propose an original and lightweight procedure in Appendix B. A principal component analysis is then performed to

estimate the basis $(h_i)$ (see Appendix A). Each PSF is projected on this basis and the resulting coefficients are interpolated spatially to provide an estimate of the functions $(\alpha_i)$, which can then be discretized as $(\alpha_i)$ (see Appendix C). All those steps are subtle and need to be performed carefully to obtain precise estimates. The technical details reported in the appendices are therefore of great importance. While revising the manuscript, we realized that the astronomical software [11], [12] was actually proposing many similar ideas for a different purpose.

As an output of the algorithm, the pair $(h_i, \alpha_i)_{1\le i\le I}$ provides a complete description of the operator, since we know an approximation of the PSFs at each image location. The integer $I$ is a user provided parameter.

### D. Estimating a subspace of operators

a) **Motivation:** A microscope produces different transfer functions depending on physical parameters that can be hard to control. Typical examples include temperature variations, focal screws, small tilts of optical elements, surface flatness of cover-slides, slight variations of a spatial light modulator rest state,... In those conditions, capturing a single operator (as proposed in the previous section) to describe the microscope might lead to model mismatches and reconstruction errors. In this section, we propose an alternative approach where we aim at learning a family of plausible operators that capture all the possible states of a microscope. The principle and the mathematical foundations behind this approach (statistical properties and fields of application) were recently established by a subset of the authors in [20]. We refer the interested reader to this paper for more details. We provide a simplified description below.

b) **Principle:** The first requirement to apply this technique is to image stacks of fluorescent micro-beads (in 2D or 3D) under multiple conditions. This process can be automatized when using advanced optical tables with motorized stages and thermostatic chambers. An alternative is to probe only the “extreme” conditions (e.g. highest and lowest plausible temperatures and tilts). After this experimental process is achieved, we have access to a set of images $(u_k)_{1\le k\le K}$. The idea of our estimation procedure is to apply the following procedure:

1) For each image $u_k$, extract the most relevant PSF patches (see Appendix A) and remove the background (see Appendix B).

2) Apply a principal component analysis to the set of patches over multiple images and keep $I$ principal components (the PSF basis).

3) Apply a z-score test to discard the patches that are likely outliers (e.g. multiple PSFs in a patch).

4) Reapply a principal component analysis to better estimate the principal components.

5) For each image $u_k$ and each coefficient $i$, interpolate the coefficient maps $\alpha_{i,k}$ (see Section C). This interpolation process is subtle: in particular we provide a
<figure><img src="image_2.png" /><figcaption>Fig. 1: Structure of the algorithm for single operator estimation. (a) Background removal procedure. (b) Selection of well isolated PSFs. (c) Extraction of relevant PSF patches. (d) Principal Component Analysis of the PSFs to find a low dimensional basis. (e) Projection of each selected PSF on the low-dimensional basis. (f) Interpolation of the PSFs coefficients using radial basis functions and correction to ensure admissible PSFs.</figcaption></figure>

novel method to learn features such as nonnegativity, or the natural decay of coefficients on the PSF basis.

6) Apply a randomized principal component analysis [31] to the whole set of sampled interpolation maps $(\alpha_{i,k})_{1 \le i \le I, 1 \le k \le K}$. It is often necessary to apply a randomized SVD <fnref n="3" /> here since the interpolation maps $\alpha_{i,k}$ are typically large images.

7) Keep the $J$ largest principal components $(a_j)_{1 \le j \le J}$.

8) Project each interpolation map $\alpha_{i,k}$ onto the basis $(a_j)_{1 \le j \le J}$, to obtain the matrices $\Gamma_k \in \mathbb{R}^{I \times J}$ defined by

$$ \Gamma_k[i, j] = \langle \alpha_{i,k}, a_j \rangle. $$

The output of this process is two orthogonal bases $(h_i)_{1 \le i \le I}$ (which describe the PSFs compactly) and $(a_j)_{1 \le j \le J}$ (which describe the PSFs variations compactly) as well as a set of matrices $(\Gamma_k)_{1 \le k \le K}$ in $\mathbb{R}^{I \times J}$ (which describe the operators associated to each image $u_k$). The

operator $H_k$ associated to the $k$-th input image $u_k$ is then defined for all $u$ by:

$$ H_k u = \sum_{1 \le i \le I} \sum_{1 \le j \le J} \Gamma_k[i, j] h_i \star (a_j \odot u). $$

c) Reducing the family of admissible operators: The subspace of operators $\mathcal{H}$ that compactly describes the possible operators is defined by

$$ \mathcal{H}^{\text{def.}} = \text{span}(\mathbf{L}_{i,j}, 1 \le i \le I, 1 \le j \le J) $$

where $\mathbf{L}_{i,j}u \stackrel{\text{def.}}{=} h_i \star (a_j \odot u)$ is a simple product-convolution operator. The dimension of $\mathcal{H}$ is $I \times J$.

However, all operators in the subspace $\mathcal{H}$ are not plausible. For instance, all PSFs are nonnegative, which is often a critical feature to avoid ringing artifacts. It is possible to further restrict the family of operators as follows. Assuming that all the extreme conditions have been explored, we can construct the *convex hull* of the coefficients $\Gamma_k$:

$$ \begin{aligned} \mathcal{C} &\stackrel{\text{def.}}{=} \text{conv}(\Gamma_k, 1 \le k \le K) \\ &\stackrel{\text{def.}}{=} \left\{ \sum_{k=1}^{K} \lambda_k \Gamma_k, \lambda_k \ge 0, \sum_{k=1}^{K} \lambda_k = 1 \right\}. \end{aligned} $$

The quality of the estimate $\mathcal{C}$ with respect to the number of observations $K$ was studied in [20]. If all the sampled

<footnote n="3">³A standard way to compute principal components requires computing a singular value decomposition. To retrieve the $k$ first components, for an $m \times n$ matrix, the complexity is $O(mnk)$, which is intractable for large scale computations. On its side, the randomized SVD provides a certified approximate solution with a complexity in $O(\log(k)mn)$ and requires significantly less memory.</footnote>
PSFs are nonnegative, then any conical combination is non-negative too and imposing the matrix $\Gamma \in C$ will therefore preserve this property. Another important feature can be preserved: the coefficients $\Gamma_k[i, j]$ follow a distribution that decays in average with $i$ and $j$, since they correspond to eigenvectors of decreasing importance. The set $C$ also captures this property, resulting in more realistic PSFs.

### E. Implementation details

a) Normalizing the operators: Of importance, let us mention that the procedure described previously suffer from a well known identifiability issue. Since the micro-beads intensity is usually unknown, the operator can be estimated only up to multiplicative constant. To address this problem, it is possible to replace the matrices $\Gamma_k$ by a normalized version

$$ \bar{\Gamma}_k = \Gamma_k / \sum_{i,j} \Gamma_k[i, j] $$

and to replace the convex hull by the conical hull

$$ C = \mathrm{cone}(\Gamma_k, 1 \le k \le K) \stackrel{\mathrm{def.}}{=} \left\{ \sum_{k=1}^{K} \lambda_k \Gamma_k, \lambda_k \ge 0 \right\}. $$

b) Normalizing the PSFs: The PSFs need to be normalized in different ways. First, they need to be registered at a subpixel accuracy. To do so, we propose a method described in Appendix A that amounts to assuming that their center of mass is located at the origin.

Second, it may also happen that the micro-beads are not perfectly identical and have different fluorescence levels. In that case it is important to normalize the PSF patches (after background removal) by imposing that they sum to 1. By doing so, the operators will be estimated without accounting for the variations of intensity that they may induce due to non uniform illumination. This effect can still be captured by assuming that the loss of intensity is proportional to the background. It then suffices to multiply the normalized patches by the background estimate.

Finally, we also normalize the $l^2$-norm of the PSFs prior to computing the principal component analysis, in order to give the same importance to every PSF.

c) Selecting the subspace sizes: The subspace sizes $I$ and $J$ are the two values that a user needs to provide in order to estimate the subspace. If the number $I$ (related to the subspace of PSFs) is too small, then the PSFs will be badly reproduced, while a value that is too large will result in noisy operators (the so-called over-fitting in machine learning). Similarly, the number $J$ captures the variations of the PSFs and has to be chosen with caution. Finally, we would like $I$ and $J$ to be as small as possible to reduce the computing times: the cost of applying a product-convolution operator is directly proportional to $I$.

The simplest way to choose $I$ is to test different values on a subset of representative PSFs and keep the lowest value that leads to a visually decent reconstruction of the PSFs. The same can be done with $J$. In practice, we observed

that the values $I = 5$ and $J = 5$ faithfully reproduce the operators from a perceptual point of view.

Another possibility is to apply recent results in statistics [32] that provide a simple and optimal way (under a Gaussian noise assumption) to choose $I$ and $J$. The rule consists in keeping the principal components associated to a singular value larger than $2.858 \cdot \sigma_{med}$, where $\sigma_{med}$ is the median of the set of singular values. This procedure requires computing the set of all singular values to evaluate the median. In practice, this is possible only for the PSF patches which are low dimensional.

### III. RESULTS

In this section, we test the proposed algorithms against 2 different data-sets: the first one is simulated while the other comes from a wide-field microscope. We start with the estimation of a single operator and of a subspace of operators. We conclude the paper with an application to the deblurring of sparse objects.

#### A. Data-sets

1) Simulation: We generate several product-convolution operators by designing a collection of admissible PSFs and space variations. The collection of PSFs is obtained by taking all the slices a 3D astigmatic PSF $h(x_1, x_2, z)$ where $z$ denotes the variations in the optical axis direction. The expression of the 2D PSF at a distance $z$ from the focal plane (used in [7] for instance) is given by:

$$ h(x, z) \propto \left| \int_D \left( \exp(-2i\pi(\langle \xi, x \rangle)E(\xi, z)) \right) d\xi \right|^2, \quad (6) $$

where $x = (x_1, x_2)$, $\xi = (\xi_1, \xi_2)$, $D$ is a disk of radius $NA/\lambda$ and $E$ is the electric field at the pupil plane given by

$$ E(\xi, z) = \exp \left( 2i\pi \left( \sum_i c_i Z_i(\xi) + z \sqrt{\frac{n^2}{\lambda^2} - |\xi|^2} \right) \right). \quad (7) $$

Here, $Z_i$ denotes the $i$-th Zernike polynomial, $NA$ is the numerical aperture, $n$ is the refractive index of the immersion oil and $\lambda$ is the emission wavelength. Note that we neglected the lateral displacements, which is a crude approximation for a large numerical apertures.

We then produce a collection of 3D PSFs by varying the parameters of the model. Each PSF is produced by taking random values of the parameters $NA \in [1.1, 1.4]$, $\lambda \in [550, 650]$, $n \in [1.41, 1.61]$. The parameters $c_i$ have been fixed to reproduce an astigmatic PSF. In addition, we model spatial variations by varying the depth $z$ in the lateral direction. We use random polynomials of low degree to generate various depth profiles. An additive background is generated with a smooth Gaussian random process, and the final image is degraded with Poisson noise, see Fig. 2a for a simulation example.
<figure><img src="image_3.png" /><figcaption>(a) Simulated micro-beads image.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Zooms on a few PSFs in Fig. 2a.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Operator estimated from the single image in Fig. 2a.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) Operator estimated from a family of 50 images.</figcaption></figure>

Fig. 2: Simulation experiment: synthesized operators are applied to randomly scattered Dirac masses in Fig. 2a (512 × 512 pixels). We then test the estimation procedure using a single micro-bead image in Fig. 2c and from 50 images in Fig. 2d. Observe that the estimation from a family is far less noisy.

<figure><img src="image_7.png" /><figcaption>(a) 2304 × 2304 image of fluorescent micro-beads</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Zooms on a few PSFs of the original image.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Estimated operator from the single image in Fig. 3a.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) Operator estimated from a family of 18 images.</figcaption></figure>

Fig. 3: Image of micro-beads taken with a wide field microscope and estimation results. The contrasts have been stretched for a better visualization. Similarly to the simulation example in Fig. 2, the operator seems far better reconstructed using the family of images.

2) *Experimental data:* In all experiments, we used a perfectly plane mono-layer of 100nm diameter micro-beads. There is no refractive index mismatch between the cover-slide and the immersion oil, allowing to avoid spherical aberrations. We used a wide-field fluorescence microscope with a ×100 objective lens (CFI SR APO 100XH NA 1,49 DT 0,12 Nikon) mounted on a Nikon Eclipse Ti-E and a Hamamatsu sCMOS camera (ORCA FLASH4.0 LT). A lens with 1.5 magnification (this is an additional magnification available on Nikon Eclipse Ti-E) is used to obtain 43nm pixel size on the image plane. A 200-nm Z interval was acquired on each image. We use SPECTRA X light engine with excitation of 633nm, and emission of 670nm. Micromanager software was used for the acquisition software. This produces images of 2304 × 2304 pixels. We collected 18 stacks of fluorescent micro-beads, each one is 8µm thick and is composed of 21 z-stacks. We keep only the 5 central slices since the beads are too degraded when far away from the focal plane. This amounts to a total of 90 images and more than 9700 2D PSFs. We display one image in Fig. 3a.

## B. Estimating operators

In this section we illustrate some features of the proposed methods by estimating a single operator and a subspace of operators from the image of micro-beads generated by the previously described microscopes. Each experiment is performed on a workstation equipped with Intel Xeon E5 and a GPU card Nvidia Tesla K20c from 2012 (2019 technologies are expected to be 4 times faster).

1) *Simulation:* We apply the proposed estimation procedures both for a single operator and for a subset of operators. The computing time for a single operator is 15 seconds when estimating 3 principal elements for the PSFs and 3 principal components to describe the coefficients variations. To estimate the subspace of simulated operators, we used 50 different micro-beads images. The computing times increased to 500 seconds (i.e. 10 second per image). The results are displayed in Fig. 4. Of importance, notice that the results obtained with the subspace of operators are based on micro-beads images generated with different operators and in particular different PSFs. Despite this higher variety of possible shapes, the method is able to automatically infer the common patterns and to achieve better denoising and estimation performance thanks to the larger number of observed PSFs.
The estimated operators can be visualized by applying them to a Dirac comb, see Fig. 2c and 2d. To compare the quality of reconstruction, we simply evaluate a rescaled $ℓ^2$ distance between the resulting images and the true one. Given two images $u_0$ and $u$, the rescaled $ℓ^2$-distance is computed as $\inf_{a \in \mathbb{R}, b \in \mathbb{R}} \|au + b - u_0\|_2$. It is compulsory to rescale the distance, since there is an inherent ambiguity between the micro-beads intensity and the microscope gain.

Estimating a subspace of operators rather than a single operator improves the quality of the reconstruction allowing to go from 50% to 15% of relative distance between the estimation and the ground truth. From a perceptual point of view, the PSFs family is no longer corrupted with noise, and the coefficients maps family seems smoother.

2) *Wide-field microscopy:* We estimate a single operator and a whole subspace of product-convolution operators based on the 2304 × 2304 images from the data-set from wide-field microscopy. We set $I = 5$ and $J = 5$. Estimating a single operator takes about 150 seconds using 120 PSFs, while estimating the whole subspace takes about 3 hours (i.e 2 minutes per operator) using 9700 PSFs. These computing times are remarkable given the computer features and given that the complete dataset weighs about 5Gb.

The PSFs and space variations bases are displayed in Fig. 5 and the estimated operators are displayed in Fig. 3. Similarly to the previous section, we observe that the basis and operators obtained using a large set of images is significantly less noisy. While the principal components beyond 3 contain a significant amount of noise for the single image, the 5-th component of the subspace approach still seems to contain useful geometrical features. The improvement of the coefficients maps is harder to evaluate since the corresponding convolution kernels have changed. Overall, learning the subspace of operators led to a significantly improved reconstruction of the operator with no visible residual noise remaining in Fig. 5. A large part of the improvement comes from the fact that more PSFs are observed and that the noise can be averaged out. The second reason is that the estimates of space variations are less sensitive to errors and turn out to be smoother.

### C. Blind deblurring

Image deblurring is a technique that can lead to significant improvements of image resolution and quality. In most acquisitions, this technique is however neglected since it requires strong skills in optics, image processing and computer science. In particular, the prior calibration of a microscope is critical: model mismatches can lead to dramatic performance losses and oftentimes lead biologists to prefer using the raw images. In this paragraph, we show that the proposed methodology of learning a whole family of operators to describe the microscope allows to avoid a precise calibration before each experiment and therefore significantly eases the application of a deblurring algorithm.

The key observation is that identifying an operator from a single degraded image becomes rather easy when the

operator depends linearly on a small number of parameters. In particular, if we know beforehand that the degraded image contains a few point sources, we show that a simple constrained least squares problem allows to recover the operator in Appendix D. We then design an original non-blind deblurring algorithm with a known operator for sparse + smooth images. The overall algorithm is called BSS for Blind Sparse + Smooth deblurring. The idea of using two steps algorithms to perform blind deblurring based on reduced models was already explored in [24], [33] for the specific case of motion deblurring. Our approach however differs significantly in the way we model the blur, in the estimation and deblurring process.

To assess the proposed methodology, we test the proposed algorithm on a real image of fluorescent micro-beads aligned along filament like structures. We perform this experiment on an image obtained with the same wide-field microscope as the one used to collect the real data. The blurry image in Fig. 6 is acquired at a distance of 400nm from the focal plane. It is possible to compare the image at the focal plane in Fig. 7, (ii) with the image used for the deblurring experiment in Fig. 7, (i).

We first compare the estimation of the PSF from the image in Fig. 6, using different approaches. In Fig. 6, (i), (ii), (iii) we show 16 equidistant PSFs estimated using 3 different approaches. In (i), we used the Matlab code `deconv.blind` based on an alternate minimization of two quadratic criteria. We set 20 iterations and initialized the method with a PSF size of 21 × 21. The PSF size is clearly underestimated. In Fig. (ii), we partitioned the image into 4 × 4 patches. Within each patch, we estimated the PSF by averaging multiple isolated PSFs as is usually recommended, see e.g. http://python-microscopy.org/doc/PSFExtraction.html. In average, we could only use two PSFs within each patch since the bead density is high and only isolated PSFs can be used. Therefore, the PSFs are still noisy, and their shape seems inaccurate, especially on the top-left corner. In Fig. (iii), we show the output of our blind identification algorithm. The PSFs are denoised and it seems that we can better reproduce the first ring of the PSF, though this effect cannot be quantified in this experiment.

We then propose comparisons with different deblurring algorithms and show the result on some patches in Fig. 7. In (i), we show the image used as an input for the deblurring. In (ii), we show the image at the focal plane. In (iii), we show the result of the function `deconv.blind` from Matlab. Here we assumed that the blur was piecewise constant on each of the 16 patches. In (iv), we show the result obtained with the software Huygens Professional version 19.04 (Scientific Volume Imaging, The Netherlands). The choice of Huygens software is motivated by its wide use among research facilities. It allows to perform a patch-wise deblurring of the full image. Here, we used the 4x4 patch decomposition in Fig. 6, (ii). The reconstruction is displayed in Fig. 7, (iii). We also conduct a second comparison with the open-source software `DeconvolutionLab2` [34]. Again,
<figure><img src="image_11.png" /><figcaption>Fig. 4: Estimation on simulated operators. We compare the estimation with a single operator (left) and with a set of 50 randomly sampled operators (right). Top: the PSF basis ($h_i$)$_{1\le i\le 3}$. Bottom: the space variations basis ($a_j$)$_{1\le j\le 3}$. Observe that the PSF basis is significantly less noisy and that the space variations are significantly smoother when estimating over the set of 50 images.</figcaption></figure>

this plugin is unable to identify the blur and we feed it with the operator estimated in Fig. 6, (ii). We use the Richardson-Lucy with total-variation regularization, which provides the most satisfactory results, see Fig. 10. The reconstruction is displayed in Fig. 7, (v). Finally, we show the output of the proposed Blind Sparse+Smooth algorithm in Fig. (vi). The Matlab `deconv.blind` approach clearly outputs unsatisfactory results (here we show the best achievable result by manual tuning of the parameters). The Huygens software was used with the default parameters. We input a 2D PSF extended in 3D to account for the defocus, specify 200 iterations with the CMLE algorithm. We also use the `DeconvolutionLab2` software with a Richardson-Lucy algorithm regularized with the total variation. Both algorithms identify single molecules, but also produce a significant amount of ringing. Finally, the output of our blind deblurring algorithm in Fig. 7, is really convincing. It rather faithfully reproduces the image obtained at the focal plane in Fig. 7, (ii) with an even better resolution. Observe that this is a really challenging setting: the input image has a significant amount of noise and while the PSFs on the left part are rather small, their diameter is about 40 pixels on the right of the field of view.

## IV. DISCUSSION

In this section, we discuss some limitations and possible extensions of the proposed approach.

a) **Centering the PSFs:** The detection of PSFs works by finding the maximum of correlation of a Gaussian with the images of micro-beads. The implicit assumption behind this procedure is that the PSFs have a center of mass located at the origin. Unfortunately, this hypothesis is wrong for some abberrations such as coma. In that case, the proposed method will result in operators that - in addition to blur - produce slight deformations of the image. Unfortunately, without prior assumption on the PSF center, it is impossible

to resolve the ambiguity between the micro-bead position and the center of mass of the PSF. In general, we can therefore expect slight distortions of the images with the proposed approach.

b) **Physicality of the PSFs:** The proposed methodology is able to reproduce some features of real point spread functions such as nonnegativity and natural decay of the coefficients in the PSF basis, thanks to the projection step on the conical hull of observed operators. This feature is important and original. When looking at the PSFs inferred by our algorithms, see e.g. Fig. 3, we can however see that they are not entirely satisfactory. For instance, the dark rings that can be seen on real diffraction patterns are not reproduced accurately. At this stage we do not know whether it is possible to obtain them with purely data driven approach as the one proposed here since they are not visible on the acquired images, which suffer from numerous artifacts such as noise, quantization, sampling, and background addition.

c) **2D versus 3D:** All the proposed algorithms and examples have been implemented in 2D, but their extension to 3D is straightforward. From a practical point of view, the estimation of a 3D operator requires to image scattered microbeads in a medium such as an agarose gel. This would allow to characterize the optical response of other types of microscopes such as confocal microscopes or light sheet fluorescence microscopes.

A limitation of the proposed approach for 2D microscopy is that the variations of the PSF in depth are important and the object is really 3D. In that case, the microscope response should be modelled as an operator mapping 3D functions to 2D images and we should infer 3D PSFs from 2D slices, which is significantly harder than what we did here. For instance, the method does not apply to 3D super-resolution microscopy.
<figure><img src="image_12.png" /><figcaption>(a) Operator estimation with a single image containing 94 micro-beads, see Fig. 3a. Top: the PSF basis ($h_i$)$_{1\le i\le 5}$. Bottom: the space variations basis ($a_j$)$_{1\le j\le 5}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Subspace estimation with 90 images and 9700 PSFs.</figcaption></figure>

<caption>Fig. 5: Learning the PSF and space variations bases for a standard wide-field microscope.</caption>

V. CONCLUSION

We proposed a set of algorithms to learn a set of product-
convolution representations of optical responses in fluores-
cence microscopy. One of the main originality is to estimate
a subspace of operators to capture the whole diversity of
possible space-varying blurs of a given microscope. This
is in sharp contrast with existing approaches which simply
characterize the microscope by a single PSF, or - at best -
by a single spatially variant operator.

An important outcome of this work is that the learned sub-
space strongly improves the identifiability in blind-inverse
problems such as blind deblurring or blind super resolution.
These arguably constitute two of the most challenging issues
in computational imaging. For instance, recent theoretical
progress based on lifting techniques [35] require the prior
knowledge of a low dimensional subspace. We have pro-
posed an original blind deblurring approach coined BSS to
efficiently solve this problem when imaging point sources
with a smooth background.

Future work will consist in extending the existing codes
to 3D images, providing an open-source toolbox together
with realistic responses of microscopes. This will enable
the optics and signal processing communities to test their
algorithms against realistic operators.

We expect the proposed work to have far reaching applica-
tions ranging from the metrology of imaging systems to new
advanced microscopy methods such as supercritical angle

localization microscopy, metal enhanced fluorescence, po-
larization microscopy. All these applications require highly
accurate models which are currently unavailable for large
fields of view.

ACKNOWLEDGMENTS

This work was supported by Fondation pour la Recherche
Médicale (FRM grant number ECO20170637521 to V.D)
and by the GDR ISIS project FiMOSuReMi. It was also
supported by the ANR JCJC Optimization on Measures
Spaces ANR-17-CE23-0013-01. The authors want to thank
Emmanuel Soubies warmly for his valuable comments.

REFERENCES

[1] M. Born and E. Wolf, *Principles of optics: electromagnetic theory of propagation, interference and diffraction of light*. Elsevier, 2013.

[2] C. Preza and J.-A. Conchello, "Depth-variant maximum-likelihood restoration for three-dimensional fluorescence microscopy," *JOSA A*, vol. 21, no. 9, pp. 1593-1601, 2004.

[3] J. W. Shaevitz and D. A. Fletcher, "Enhanced three-dimensional deconvolution microscopy using a measured depth-varying point-spread function," *JOSA A*, vol. 24, no. 9, pp. 2622-2627, 2007.

[4] G. Zheng, X. Ou, R. Horstmeyer, and C. Yang, "Characterization of spatially varying aberrations for wide field-of-view microscopy," *Optics express*, vol. 21, no. 13, pp. 15131-15143, 2013.

[5] A. von Diezmann, M. Y. Lee, M. D. Lew, and W. Moerner, "Correcting field-dependent aberrations with nanoscale accuracy in three-dimensional single-molecule localization microscopy," *Optica*, vol. 2, no. 11, pp. 985-993, 2015.

[6] J. W. Goodman, *Introduction to Fourier optics*. Roberts and Company Publishers, 2005.
<figure><img src="image_14.png" /><figcaption>Fig. 6: Large 2048 × 2048 crop of the original image taken at 400nm from the focal plane. The contrasts have been stretched for a better visualization. In (i), (ii) and (iii), we show the PSFs estimated using different approaches at the centers of the 4 × 4 uniform tiling of the image represented in gray. In (i), we used the *deconv blind* algorithm from Matlab. In (ii), we averaged isolated PSFs within each patch. In (iii), we show the output of our blind identification algorithm.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 7: Zooms on the output of different deblurring algorithms. (i) image used for the deblurring experiment. (ii) image at the focal plane. (iii) output of the Matlab *deconv blind* algorithm. (iv) output of the Huygens software. (v) output of *DeconvolutionLab2* with a Richardson-Lucy TV algorithm. (vi) output of our Blind Sparse+Smooth algorithm.</figcaption></figure>

[7] A. Aristov, B. Lelandais, E. Rensen, and C. Zimmer, "Zola-3d allows flexible 3d localization microscopy over an adjustable axial range," *Nature communications*, vol. 9, no. 1, p. 2409, 2018.

[8] T. Yan, C. J. Richardson, M. Zhang, and A. Gahlmann, "Computational correction of spatially variant optical aberrations in 3d single-molecule localization microscopy," *Opt. Express*, vol. 27, no. 9, pp. 12582–12599, Apr 2019.

[9] S. F. Gibson and F. Lanni, "Diffraction by a circular aperture as a model for three-dimensional optical microscopy," *JOSA A*, vol. 6, no. 9, pp. 1357–1367, 1989.

[10] J. Anderson and I. R. King, "Toward high-precision astrometry with wfpc2. i. deriving an accurate point-spread function," *Publications of the Astronomical Society of the Pacific*, vol. 112, no. 776, p. 1360, 2000.

[11] E. Bertin and S. Arnouts, "Sextractor: Software for source extraction," *Astronomy and Astrophysics Supplement Series*, vol. 117, no. 2, pp. 393–404, 1996.

[12] E. Bertin, "Psfex: Point spread function extractor," *ascl*, pp. ascl–1301, 2013.

[13] Y. Li, M. Mund, P. Hoess, J. Deschamps, U. Matti, B. Nijmeijer, V. J. Sabinina, J. Ellenberg, I. Schoen, and J. Ries, "Real-time 3d single-

molecule localization using experimental point spread functions," *Nature methods*, vol. 15, no. 5, p. 367, 2018.

[14] M. Elad and Y. Hel-Or, "A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur," *IEEE Transactions on Image Processing*, vol. 10, no. 8, pp. 1187–1193, 2001.

[15] F. N. Mboula, J.-L. Starck, S. Ronayette, K. Okumura, and J. Amiaux, "Super-resolution method using sparse regularization for point-spread function recovery," *Astronomy & Astrophysics*, vol. 575, p. A86, 2015.

[16] M. Gentile, F. Courbin, and G. Meylan, "Interpolating point spread function anisotropy," *Astronomy & Astrophysics*, vol. 549, p. A1, 2013.

[17] C. Chang, P. Marshall, J. Jernigan, J. Peterson, S. Kahn, S. Gull, Y. AlSayyad, Z. Ahmad, J. Bankert, D. Bard et al., "Atmospheric point spread function interpolation for weak lensing in short exposure imaging data," *Monthly Notices of the Royal Astronomical Society*, vol. 427, no. 3, pp. 2572–2587, 2012.

[18] J. Bigot, P. Escande, and P. Weiss, "Estimation of linear operators from scattered impulse responses," *Applied and Computational Harmonic Analysis*, vol. 47, no. 3, pp. 730–758, 2019.

[19] F. Ngolè, J.-L. Starck, K. Okumura, J. Amiaux, and P. Hudelot, "Constraint matrix factorization for space variant psfs field restoration," *Inverse Problems*, vol. 32, no. 12, p. 124001, 2016.

[20] V. Debarnot, P. Escande, and P. Weiss, "A scalable estimator of sets
of integral operators,” Inverse Problems, 2019.

[21] A. Pinkus, N-widths in Approximation Theory. Springer Science & Business Media, 1985, vol. 7.

[22] F. Heide, M. Rouf, M. B. Hullin, B. Labitzke, W. Heidrich, and A. Kolb, “High-quality computational imaging through simple lenses,” ACM Transactions on Graphics (TOG), vol. 32, no. 5, pp. 1–14, 2013.

[23] R. C. Flicker and F. J. Rigaut, “Anisoplanatic deconvolution of adaptive optics images,” J. Opt. Soc. Am. A, vol. 22, no. 3, pp. 504–513, Mar 2005.

[24] M. Hirsch, C. J. Schuler, S. Harmeling, and B. Schölkopf, “Fast removal of non-uniform camera shake,” in 2011 International Conference on Computer Vision. IEEE, 2011, pp. 463–470.

[25] P. Escande and P. Weiss, “Approximation of integral operators using product-convolution expansions,” Journal of Mathematical Imaging and Vision, vol. 58, no. 3, pp. 333–348, 2017.

[26] L. Denis, E. Thiébaut, F. Soulez, J.-M. Becker, and R. Mourya, “Fast approximations of shift-variant blur,” International Journal of Computer Vision, vol. 115, no. 3, pp. 253–278, 2015.

[27] M. Arigovindan, J. Shaevitz, J. McGowan, J. W. Sedat, and D. A. Agard, “A parallel product-convolution approach for representing depth varying point spread functions in 3d widefield microscopy based on principal component analysis,” Opt. Express, vol. 18, no. 7, pp. 6461–6476, Mar 2010.

[28] N. Patwary and C. Preza, “Image restoration for three-dimensional fluorescence microscopy using an orthonormal basis for efficient representation of depth-variant point-spread functions,” Biomedical optics express, vol. 6, no. 10, pp. 3826–3841, 2015.

[29] F. Sroubek, J. Kamenicky, and Y. M. Lu, “Decomposition of space-variant blur in image deconvolution,” IEEE Signal Processing Letters, vol. 23, no. 3, pp. 346–350, 2016.

[30] D. Miraut and J. Portilla, “Efficient shift-variant image restoration using deformable filtering (part i),” EURASIP Journal on Advances in Signal Processing, vol. 2012, no. 1, p. 100, 2012.

[31] N. Halko, P.-G. Martinsson, and J. A. Tropp, “Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions,” SIAM review, vol. 53, no. 2, pp. 217–288, 2011.

[32] M. Gavish and D. L. Donoho, “The optimal hard threshold for singular values is $4/\sqrt{3}$,” IEEE Transactions on Information Theory, vol. 60, no. 8, pp. 5040–5053, 2014.

[33] O. Whyte, J. Sivic, A. Zisserman, and J. Ponce, “Non-uniform deblurring for shaken images,” International journal of computer vision, vol. 98, no. 2, pp. 168–186, 2012.

[34] D. Sage, L. Donati, F. Soulez, D. Fortun, G. Schmit, A. Seitz, R. Guiet, C. Vonesch, and M. Unser, “DeconvolutionLab2: An open-source software for deconvolution microscopy,” Methods—Image Processing for Biologists, vol. 115, pp. 28–41, February 15, 2017.

[35] A. Ahmed, B. Recht, and J. Romberg, “Blind deconvolution using convex programming,” IEEE Transactions on Information Theory, vol. 60, no. 3, pp. 1711–1732, 2013.

[36] M. D. Buhmann, Radial basis functions: theory and implementations. Cambridge university press, 2003, vol. 12.

[37] Y. Nesterov, Lectures on convex optimization. Springer, 2018, vol. 137.

[38] A. Kumar, V. Sindhwani, and P. Kambadur, “Fast conical hull algorithms for near-separable non-negative matrix factorization,” in International Conference on Machine Learning, 2013, pp. 231–239.

[39] R. Henrion and A. Seeger, “On properties of different notions of centers for convex cones,” Set-Valued and Variational Analysis, vol. 18, no. 2, pp. 205–231, 2010.

# APPENDIX A
## PSFS SELECTION AND PROCESSING

a) **Rough PSF detection:** The rough detection of PSFs on an image $u$ is quite easy: we simply identify the maxima of the Laplacian of a blurred version of $u$ and keep those above a user-provided threshold.

b) **Rejection of bad patches:** The analytical expression of a PSF yields functions that are not compactly supported. However, they decay quite fast at infinity (an Airy pattern decays as $1/|x|^3$) and the information brought by the image

is dominated by noise far away from the PSF center. In what follows, we let $r$ denote an upper-bound on the maximal radius of a PSF. To avoid patches containing more than one PSF or PSFs too close to the image boundaries, we select the subset of PSFs with a center being at least $r$ pixels from the image boundary and $2r$ pixels from the other centers. At the end of this procedure, we obtain a set of patches containing isolated PSFs.

c) **Shift and re-sampling:** The maxima of the discrete Laplacian are localized on the pixel grid and might not correspond to the PSF center. In addition, there is no reason for the micro-beads to be perfectly centered on a pixel. To obtain a better PSF centering, we find the maxima of correlation with a continuous Gaussian function. This allows us to re-interpolate the patches around this center using a bi-cubic interpolation. The output is a sequence of patches $(p_m)_{1 \le j \le M}$.

d) **Compact supports:** To avoid learning noise, we select two radii $r_{ext} \ge r$ beyond which the PSF image is dominated by noise. In practice, we set $r_{ext} = 1.2r$. We then apply a multiplier to each patch that smoothly decays to 0 from $r$ to $r_{ext}$ (a bump function). The analytical expression of the radial multiplier function $m(t)$ is given by

$$m(t) = \begin{cases} 1 & \text{if } t < r, \\ 0 & \text{if } t > r_{ext}, \\ \exp\left(\frac{-1}{1-\left(\frac{t-r}{r_{ext}-r}\right)^2} + 1\right) & \text{otherwise.} \end{cases}$$

e) **Principal component analysis:** Once all selected PSFs have been treated, we perform a principal component analysis to obtain an optimal representation basis. Depending on the number of sampled PSFs, we can use a standard singular value decomposition or a randomized one [31].

f) **Finer rejection of false detections:** Finally, it may happen that the initial selection of patches did not screen out all unwanted PSFs (e.g. too degraded or multiple ones on a patch). To get a finer selection and avoid biases, we project all the PSFs on the learned basis $(h_i)_{1 \le i \le I}$. We then remove the coefficients in $\mathbb{R}^I$ obtained by the projection which are likely outliers. We simply use a z-score test: we evaluate the empirical mean $\mu_i$ and variance $\sigma_i$ of each coefficient and discard the coefficients that deviate by more than $3\sigma_i$ from their mean. Once the unlikely patches have been removed we recompute a PCA on the remaining patches to get a better estimate of $(h_i)_{1 \le i \le I}$.

# APPENDIX B
## BACKGROUND REMOVAL

Getting a correct estimate of the background is a critical step since not accounting for it would strongly bias the PSFs estimates. The background is due to the auto-fluorescence of the sample and may vary spatially due to a non homogeneous illumination. In general this variation is smooth, allowing for an efficient estimation described hereafter.

We assume that the PSF is dominated by noise in a domain $\omega$ made of pixels outside a disk of radius $r$ from
the patch center (see the previous section for a definition of the patches). In this region, the image signal is therefore constituted of background and noise only. Assuming that the background is a smooth component, we can try to fit a low degree polynomial to these pixels. Letting $p_j$ denote the j-th patch, this amounts to solving the following quadratic problem:

$$ \inf_{\lambda \in \mathbb{R}^P} \frac{1}{2} \|\mathbf{p}_j - \mathbf{M}\lambda\|_{l^2(\omega)}^2, \quad (8) $$

where $\mathbf{M} = [\mathbf{m}_1, \dots, \mathbf{m}_p] \in \mathbb{R}^{|\omega| \times P}$ is a matrix containing the sampled low degree monomials $\mathbf{m}_p$ and $\lambda$ represents the coefficients of the polynomial. This problem boils down to a low dimensional linear system which can be solved with a linear conjugate gradient algorithm. Letting $\lambda^*$ denote the solution of this problem, the estimated background is simply $\mathbf{M}\lambda^*$. In practice, we simply use polynomials of order 2. In 2D, this yields the value $P = 6$ for the monomials 1, $x$, $y$, $xy$, $x^2$, $y^2$.

## APPENDIX C
### ESTIMATING SPACE VARIATIONS

a) **Thin-plate approximation:** Once a basis $(\mathbf{h}_i)_{1 \le i \le I}$ is computed (see section A), it is possible to project each noisy patch on this basis to get a low dimensional representation of the selected PSFs. This provides an estimate $\beta_{i,m} = \langle \mathbf{p}_m, \mathbf{h}_i \rangle$ of the values $\alpha_i(\mathbf{z}_m)$. In order to estimate the space variations, we can use surface fitting techniques on the set $(\mathbf{z}_m, \beta_{i,m})_{1 \le m \le M}$ to get an approximation of the functions $\alpha_i$.

There exist numerous surface fitting techniques. Following the numerical experiment conducted in [16], it seems that the use of radial basis function [36] is significantly more efficient than other approaches in the context of astronomy. We therefore resort to this technique.

Radial basis functions approximation can be interpreted as a variational problem in the framework of Reproducible Kernel Hilbert Spaces. In this context, the estimators $\hat{\alpha}_i$ of $\alpha_i$ can be expressed as

$$ \hat{\alpha}_i = \underset{\alpha \in H^2(\mathbb{R}^2)}{\operatorname{argmin}} \frac{1}{2} \sum_{m=1}^{M} w_m |\alpha(\mathbf{z}_m) - \beta_{i,m}|^2 + \frac{\eta}{2} |\alpha|_{H^2}^2, \quad (9) $$

where $|\alpha|_{H^2} \stackrel{\text{def.}}{=} \langle \Delta u, \Delta u \rangle_{L^2(\mathbb{R}^d)}$ and where $\eta > 0$ is a parameter that allows to trade off the proximity to the samples $\beta_{i,m}$ for the smoothness of the surface. In order to balance the importance of each PSF in the approximation, the weights $w_m$ are chosen equal to the area of the Voronoi cell associated to each location $\mathbf{z}_m$.

The solution of (9) is known to be a thin-plate spline [21] and can be computed by solving a $(M+3) \times (M+3)$ linear system.

In what follows, we will let $\alpha_i$ or $\hat{\alpha}_i$ denote a version of $\alpha_i$ sampled on a Euclidean grid.

b) **Enforcing realistic PSFs:** There is no reason for the thin-plate approximation method to generate realistic PSFs everywhere in the field of view. Indeed, the coefficients are interpolated *independently* of each other while there exists strong dependencies between them. In practice we observed that the previous method was not good at extrapolating the PSFs outside of the convex hull of the sampled PSFs. Important features like positivity for instance might be lost far away from the sampled PSFs. To avoid this effect, we propose an original framework below.

The procedure for estimating the PSF basis $(\mathbf{h}_i)$ also yields the projected coefficients $(\beta_{i,m})_{1 \le i \le I, 1 \le m \le M}$. We propose to define the set of *admissible* PSFs coefficients as

$$ \mathcal{B} \stackrel{\text{def.}}{=} \text{cone}(\beta_{\cdot,m}, 1 \le m \le M) \subset \mathbb{R}^I. $$

This roughly amounts to say that admissible PSFs correspond to the conical hull of the already observed and denoised PSFs. Taking the conical hull seems natural: if a PSF is in the set, all its scaled versions by a non-negative factor also belong to the set. Let $\alpha = (\alpha_1, \dots, \alpha_I) \in \mathbb{R}^{I \times N}$ denote a possibly infeasible estimate of interpolation map. We can generate a feasible one $\hat{\alpha}$ by projection: $\hat{\alpha}[n] \stackrel{\text{def.}}{=} \Pi_B(\alpha[n])$ for all $1 \le n \le N$.

c) **The projection algorithm:** Let $\mathcal{B} = [\beta_{\cdot,1}, \dots, \beta_{\cdot,M}] \in \mathbb{R}^{I \times M}$ denote the matrix of observed coefficients. Projecting the coefficients $\alpha \in \mathbb{R}^I$ of a PSF onto $\mathcal{B}$ amounts to solving the following convex variational problem:

$$ \inf_{\lambda \in \mathbb{R}^M, \lambda \ge 0} \frac{1}{2} \|\mathcal{B}\lambda - \alpha\|_2^2. $$

It can be solved using an accelerated projected gradient descent [37].

Unfortunately, if $M$ is very large, applying matrix-vector products with $\mathcal{B}$ for every pixel $n \le N$ becomes untractable and we need to simplify the cone $C$. Following [38], we propose to select a small subset of the columns of $\mathcal{B}$ using a simple greedy algorithm. We start with a matrix $\hat{\mathcal{B}}$ containing a single vector equal to the average of the columns of $\mathcal{B}$. We then update it by iteratively adding the column in $\mathcal{B}$ which maximizes the angle with the current conical hull of the columns in $\hat{\mathcal{B}}$. We stop when the angle is below a given threshold. In our experiments with $M = 14000$ PSFs and $I = 5$, we could obtain a very good approximation of the hull with only 20 components instead of 14000, making the projection algorithm relevant even for very large images.

## APPENDIX D
### BSS-DEBLURRING

Here, we propose a method called BSS-deblurring, where BSS stands for Blind Sparse+Smooth. Given a blurred image $u_0$, the method provides an estimate of the associated operator and a deblurred image. It consists of two separate steps: first the operator is estimated using isolated micro-beads in the image. This estimate is then used inside an original variational problem.
a) Identifying the operator: We assume that the user is able to select P patches out of $u_0$ supported on $(\omega_p)_{1\le p\le P}$ that contain isolated micro-beads. The patches $p_p$ are assumed to be well centered and without background, which can be achieved using the previously proposed methods.

Our aim is to estimate a discrete operator $H: \mathbb{R}^N \to \mathbb{R}^N$ of the form:

$$H u = \sum_{1 \le i \le I} \sum_{1 \le j \le J} \Gamma[i, j] h_i \star (a_j \odot u)$$

where the pairs of orthogonal bases $((h_i)_i, (a_j)_j)$ are known and the coefficients $(\Gamma[i, j])_{i,j}$ are unknown.

By projecting the patch $p_p$ onto the PSFs basis $(h_i)_i$, we obtain the coefficient $c_{p,i} = \langle p_p, h_i \rangle$ which is a noisy estimate of the interpolation map $\alpha_i$ at the position $z_p$:

$$H \delta_{z_p} = \sum_{1 \le i \le I} \sum_{1 \le j \le J} \Gamma[i, j] \alpha_i(z_p) h_i \approx \sum_{1 \le i \le I} c_{p,i} h_i.$$

We propose to identify $\Gamma$ by solving the following bi-linear inverse problem:

$$\underset{\substack{\Gamma \in \mathcal{C} \\ \langle \Gamma, \Gamma_0 \rangle = 1 \\ g \in \mathbb{R}_+^P}}{\operatorname{argmin}} \frac{1}{2} \sum_{1 \le p \le P} \left\| \sum_{1 \le i \le I} \sum_{1 \le j \le J} g_p \Gamma_{i,j} \alpha_j - c_{p,i} \right\|_2^2, \quad (10)$$

where $g_p$ is the unknown amplitude of the bead at position $z_p$ and $\mathcal{C}$ is the conical hull of the sampled operators. The additional linear constraint $\langle \Gamma, \Gamma_0 \rangle = 1$ is related to an intrinsic identifiability problem in blind deblurring problems: the operator can be multiplied by a constant factor and the signal by its inverse, leading to the same result. Letting $\Gamma_0$ denote a reference vector in $\mathcal{C}$, we can avoid this caveat. A nice geometrical choice consists in choosing the so-called *circumcenter* of the cone [39]. We do not discuss this choice further since the proposed method is in essence heuristic.

We solve this problem using an alternating minimization algorithm: we first solve the problem w.r.t. $g$ with fixed $\Gamma$ and then solve the problem w.r.t. $\Gamma$ with fixed $g$. The individual minimization problems are convex and can be solved with accelerated projected gradient descents.

If the amplitudes $g$ were known, the problem (10) would boil down to a constrained least square problem of size $I \times J$. Since each patch $p_p$ yields $I$ coefficients, the condition $P \ge J$ should be enough to ensure the identification. It is remarkable that such a low value (typically 5) is enough to identify the operator! Higher values of $P$ would however make the method more robust to noise.

To validate the proposed approach, we randomly select an operator in the conical hull $\mathcal{C}$, and apply it to a grid of 25 Dirac masses (for a field of view of 2304 x 2304 pixels). A significant amount of noise is added to the image and the true locations of the beads are randomly perturbed (with Gaussian random variable of variance 0.5). We then estimate the operator and display the result in Fig. 8.

We finally apply the method to the image used in the blind-deblurring example in Fig. 6. We display the reconstructed operator in Fig. 9.

<figure><img src="image_16.png" /><figcaption>Fig. 8: Identification of an operator in simulated images. a) Noisy crops used for identification. b) Impulse responses of the true operator at some locations. c) Estimated operator. Notice that the method is able to denoise the PSFs very efficiently.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 9: Identification of a real operator with 38 micro-beads. Top: 3 patches used for the estimation. Bottom: estimated PSFs.</figcaption></figure>

b) Sparse+Smooth-deblurring: Let $u_0 \in \mathbb{R}^N$ denote a blurry image and $H: \mathbb{R}^N \to \mathbb{R}^N$ a discrete integral operator. We aim at deblurring an image composed of a sparse component $u_1$ (e.g. scattered micro-beads), and a smooth component $u_2$ (e.g. auto-fluorescent background). In order to recover $u_1$ and $u_2$, we propose to solve the following original variational problem:

$$\inf_{u_1 \in \mathbb{R}_+^N, u_2 \in \mathbb{R}_+^N} \frac{1}{2} \|H(u_1 + u_2) - y\|_2^2 + \gamma_1 \|u_1\|_1 + \gamma_2 \|\Delta u_2\|_2^2.$$

The term $\frac{1}{2}\|H(u_1 + u_2) - y\|_2^2$ is the data fitting term, the term $\gamma_1 \|u_1\|_1$ promotes the sparsity of the $u_1$ component
and the term $\gamma_2 \|\Delta u_2\|_2^2$ promotes the smoothness of $u_2$. The non-negative parameters $\gamma_1$ and $\gamma_2$ allow to balance each term and have been tuned manually so as to obtain a visually pleasant result. This problem can be solved efficiently using accelerated proximal gradient descent algorithm [37].

To test the proposed algorithm, we evaluate its performance on synthetic micro-beads images with a spatially invariant Gaussian PSF with variance $\sigma = 10^{-2}$. We add a second order polynomial to simulate the background and random white Gaussian noise. The blurry-noisy image is displayed in Fig. 10a. The value of the proposed methodology comes from two distinct features: a more accurate model of microscope and a better deblurring model with the Sparse+Smooth prior. To disentangle the respective contributions of each aspect, we conduct two experiments.

We first show the impact of an accurate model. We apply the Sparse+Smooth algorithm with a PSF smaller than the true one ($\sigma = 0.5 \times 10^{-2}$) in Fig. 10d, with the true PSF ($\sigma = 10^{-2}$) in Fig. 10e, and with a PSF larger than the true one ($\sigma = 2 \times 10^{-2}$) in Fig. 10d. As can be seen in Fig. 10e, 10f and 10d, only the algorithm run with the correct PSF is able to correctly localize all sources, even when they are close together. This shows the importance of describing the microscope response accurately.

Second, we compare the Sparse+Smooth deblurring algorithm with other standard approaches implemented in DeconvolutionLab2 [34] using the true PSF ($\sigma = 10^{-2}$). We have selected the following popular methods: the regularized inverse filtering in Fig. 10g, the Richardson-Lucy TV algorithm in Fig. 10h, and FISTA algorithm with $l^1$ penalization of the Haar wavelet coefficients in Fig. 10g. In all experiments we tuned the parameters manually so as to obtain the best results from a perceptual point of view. The Sparse+Smooth algorithm seems to be by far the preferable approach to detect point sources over a smooth background.

<figure><img src="image_18.png" /><figcaption>Fig. 10: Validation of the Sparse+Smooth deblurring algorithm on a synthesized image. We first challenge the algorithm with inaccurate models of PSFs (the exact PSF, a smaller one and a larger one). As can be seen, only the true PSF produces near exact results, but the algorithm behaves nicely even with some inaccuracies. We also compare the Sparse+Smooth deblurring algorithm with three popular methods proposed in DeconvolutionLab2 [34]. Here, the proposed model performs significantly better.</figcaption></figure>