Solving the boundary artifact for the enhanced deconvolution algorithm SUPPOSe applied to fluorescence microscopy

The SUPPOSe enhanced deconvolution algorithm relies in assuming that the image source can be described by an incoherent superposition of virtual point sources of equal intensity and finding the number and position of such virtual sources. In this work we describe the recent advances in the implementation of the method to gain resolution and remove artifacts due to the presence of fluorescent molecules close enough to the image frame boundary. The method was modified removing the invariant used before given by the product of the flux of the virtual sources times the number of virtual sources, and replacing it by a new invariant given by the total flux within the frame, thus allowing the location of virtual sources outside the frame but contributing to the signal inside the frame.


Introduction
The challenge of obtaining high-resolution images beyond the limit imposed by diffraction (resolution beyond the Abbe limit), aberrations and noise with a single image was pursued for decades by deconvolution of the image with the instrument response function or point spread function (PSF). Recently in [1] the authors examine several standard deconvolution algorithms: Regularized inverse filter, Tikhonov regularization, Landweber, Tikhonov-Miller, Richardson-Lucy, and fast iterative shrinkage-thresholding. These algorithms are evaluated for different 3D-images and implemented on an opensource software. For all these deconvolution methods a marginal resolution improvement was obtained and in some cases the retrieval has shown different type of artifacts. The reason for the lack of significant resolution improvement and the presence of artifacts is that the deconvolution problem, i.e finding the function R (x, y) that convolved with the PSF yields the measured image, is that due to the noise in the measurement the deconvolution is ill posed. This happens when we use the naive method of minimizing the l 2 norm, which result in spurious oscillations or as in iterative Landweber method [2] which tends to produce an over-fitting on the noise when the number of iteration is not chosen appropriately (the number of iterations is a pseudo regularization parameter) or in Richardson-Lucy [3,4] that typically produces ringing artifacts (amplification of the noise). In fact, after deconvolution, solutions are found with negative values for R even when we are dealing with an intensity signal that must be positive. Further, Tikhonov [5], regularized inverse filtering [6] and Tikhonov-Miller methods incorporate a regularization term to convert the problem into a well posed one. There are also fast iterative methods that incorporates a sparsity constrains in the wavelet domain and where the optimization problem is solved using new parameters (step side) and a soft-thresholding operator [7]. All these regularization terms always end up being a damping factor that smooths the recovered function R. Finally, the Richardson-Lucy method with total-variation regularization method [8] add a regularization term using the l 1 norm, instead of l 2 used in the regularized inverse filtering method, to preserve image discontinuities. This l 1 term might cause artifacts if the regularization parameter is not chosen adequately producing a bad interpretation of the biological structure. These schemes have proven useful in 3D reconstruction where out of plane fluorescence can be removed but as mentioned, do not lead to a significant improvement in resolution.
On the other hand, in the deconvolving compressed sensing methods [9,10] the sparsity prior is a very restrictive constraint as shown in [11]. The authors show a good quality of the reconstruction only if there are less than three sources overlapping within the PSF. Further, in that case the authors recover R assuming that there is no background added in the image.
In the work [12] the authors define a new enhanced deconvolution method from a single image that could be applied for fluorescent 2D microscopy images even when there is an unknown background. This background includes out of focus or surrounding media fluorescence, dark counts (average), readout offset or any other additive contribution to the explored signal. The method called SUPPOSe (from superposition of virtual point sources) consists of introducing the ansatz that the objective function R (x, y) can be approximated by a sum of N virtual point sources, each one providing equal flux . The method has been tested on synthetic images showing a high increase in the resolution. In [13], the SUPPOSe method was applied to recover with high resolution cellular structures such as microtubules and F-actin filaments observed with a standard microscope with LED illumination. In these examples details down to 50 were retrieved. The ansatz introduced converts the inverse problem of finding the function R to a direct problem of finding the position of the sources, thus removing the negative valued solutions for R and as proved in [12] converting the mathematical problem to a well posed one with a resolution that depends on the number of the virtual sources selected, the quality of the PSF measurement, and mainly the signal to noise ratio. Although the resolution depends on how dense the structure is, the SUPPOSe method allows solving denser structures than in [11]. In fact, the number of optimal virtual sources and the spatial resolution achieved can be estimated and depend on the density of the structure (see supplementary material, [12]), so we do not have to add a constrain to the minimization problem to impose the sparsity condition of the image as in the compress sensing method. Since the method better suited for filiform structures such as frequently found in biological problems, and less adequate for dense volume conglomerate of emitters [12] in [14] a modified version detecting the object edges was developed.
One of the advantages of SUPPOSe is that we only have to find the positions of these virtual sources that realize a minimum of a selected fitting function. Because a background would be a dense signal an image with background is not adequate for methods assuming sparsity of the solution. To solve this problem in [12] we needed to redefine the fitting function in the cases where there was an unknown background by subtracting the average value of the image and fitting with a PSF with null average. To find the minimum of this new fitting function we have chosen a parallelizable algorithm and we observed that one way to exploit this property was tiling the image. When applying SUPPOSe to some subfigures of synthesized images we observe the presence of artifacts near the boundary. This artifact arises from the presence of fluorescent molecules close enough to the boundary (form inside and from outside the tile). This was not taken into account in the original version of SUPPPOSe [12] where  N was assumed equal to the sum of the signal within the frame. To remove these artifacts in this paper we modified the optimization problem, changing the fitting function by a new one. As in [12] we first test our method for a synthetic image, and we show that we can remove the boundary artifact present on the previous method.
Further we observe how the artifacts disappear when tiling the images.
It is important to mention that a new field has been evolving in the last years known in the literature as superresolution or high resolution image reconstruction [15]. Some methods are based in scanning beams using reversible saturable optical fluorescence transitions (RESOLFT) microscopy either by stimulated emission depletion (STED) [16,17] microscopy, ground state depletion (GSD) microscopy [18,19]. Another alternative using wide field imaging are structured light illumination (SIM) [20] that requires the acquisition of at least nine frames for an increase in resolution of a factor of 2 and that can be increased when saturation of the excited state is reached [21 -24] at the expense of even more frames. Another approach using Fresnel incoherent correlation holography (FINCH) microscopy by generating self interference using a multifocal lens [25,26] achieves a 2 to 3 times improvement and requires at least 3 images per frame. It is worthwhile pointing out that all these methods could be used in conjunction with the SUPPOSe algorithm to increase the resolution even further.
At this point a distinction must be made with the super-resolution methods for fluorescent microscopy based in the localization of single fluorescent molecules such as PALM [27,28] or STORM [29,30]. These types of methods are known as single-molecule localization microscopy (SMLM) and consist in the localization of real sources as opposed to virtual sources used in SUPPOSe. In those cases, it is basically required that the molecules do not overlap within the PSF of the microscope for a correct localization. Therefore these methods (based in high sensitivity camera detection) need hundreds or even thousands of frames to complete a single image. These methods are evolving fast and providing new insight in biological imaging [31 -34] with great effort placed in developing adequate fluorophores [35]. Much progress has been made in recent years in evolving to 3D resolution by adding to the lateral localization characteristic of STORM methods a specific technique for axial resolution. A variety of schemes have been developed including introducing an axial dependence of the PSF [36,37], using sophisticated interferometric techniques [38,39], self interferometric technique [40], supercritical-angle fluorescence [41] and more sophisticated near field schemes such as metal induced energy transfer [42]. To increase the speed methods resolving several simultaneously emitting molecules within the PSF have been implemented. This has been achieved deconvolving each image using compressed sensing so the sparsity prior is a very restrictive constraint, and this has permitted a moderate increase in STORM reconstruction speed [43 -46]. Other schemes have been approached with similar results [47,48]. As mentioned before the SUPPOSe algorithm should not be seen as a localization reconstruction but as a deconvolution operation that will not provide an extraordinary reso-lution of the SMLM methods but has the advantage of working with a single image.
The paper is organized as follows. In Section "SUP-POSe method and modified SUPPOSe method" we describe the two methods, how the fit function was adapted and which steps in the algorithm need to be modified. In Section "Results" we apply the modified SUPPOSe method to synthesized images and to real images. First we reconstruct two synthesized images using tiles by the original SUPPOSe and by the modified SUPPOSe method. In both synthetic images we observe that the modified SUPPOSe method removes the artifacts that appear in the intersection of the tiles when applying the original method. These artifacts are removed maintaining the same resolution that we had obtained when applying SUPPOSe to images supported inside the image frame. Then we use our method to remove artifacts for real images obtained from a standard fluorescence microscope. To test the performance of both methods on real images we acquired an image of mitochondrial networks using an objective with low numerical aperture (NA) and compared both reconstruction with the same sample captured with an objective with larger numerical aperture. We split the low NA image into two tiles, and we apply both methods to show that also here the modified SUPPOSe method removes the artifacts present in the old version. As well, we find that the resolution is at least that of the high NA objective. Finally we show that with the modified method we can remove artifacts of tiled F-actin wide-field images taken with high NA objectives. In this case we could distinguish structures that are distorted when applying the standard SUPPOSe method.

SUPPOSe method and modified SUPPOSe method
To explain changes introduced in the fitting function we are going to start with a brief description of the SUP-POSe method. All image measurements (in dimension D) are blurred and distorted by the Point Spread Function (PSF). We will denote the target function by R (x), the measured data S (x) and the PSF by I, with x  R D . The relation between the target function and the measured data are given by a convolution with the PSF plus the background and the noise i.e: The signal R*I and the background B are the average value of the associated stochastic variable, and the fluctuations around such averages obtained in each measurement is the noise . Hence the average value for the noise defined in this manner is null. For the case of fluorescent microscopy the background includes out of focus or surrounding media fluorescence, dark counts (average), readout offset or any other additive contribution to the explored signal. The function I is obtained after pixelation of the original PSF. Also S is only sampled for certain values 1 { } n i i x  , which are the pixels, being each x i a vector in R D where n is the number of pixels with side distance between neighbor pixels equal to d p .
The SUPPOSe method consists in approximating the target function R, by a superposition of virtual point sources of identical fluxes  so that the only unknown are the positions of the sources. That is, the approximate so- Here for each k = 1,, N, are the positions of the virtual sources, where N is the number of virtual point sources used for the fit. Observe that these positions can be repeated.
The quality of the reconstruction depends on the number of virtual sources N used. The accuracy in the positions of the virtual sources was defined in [12] by . For the precise definition of  see [12]. The authors also defined N op as the number of virtual sources N that minimizes and the corresponding optimal  was denoted by  op .

Case with no background -SUPPOSe
In the SUPPOSe method when there was no background the recorded signal S could be reconstructed approximately by Here I  was some approximation of the PSF function I (see Section 2 of the supplementary material in [12]). Given N and , we search for the position of the virtual point sources 1 { } N k k a  , that yield a minimum of the fitting function, where |||| is the standard 2-norm. Remember that D was the dimension of the space and N the number of virtual point sources used for the fit. In this case we were using that I was normalized and that the positions of the virtual point sources k a  were far from the boundary. We had in that case that, We now give a summary of the steps of the algorithm, for more details see [12]. Algorithm 1. To process an image, 1. Start with some arbitrary N and use (5).
2. Then use a Genetic Algorithm to minimize  2 in (4). 3. Make an histogram of the solution. Compute N op and  is scaled accordingly. Return to step (2). 4. Convolve the obtained virtual point sources with the known shape of the source used for the deter-mination of the PSF using  op as its width. This step is optional. The third step it is necessary when it is not easy to calculate the optimal N.

Case with background -SUPPOSe
In the case where there was an unknown background level in (1) we defined Finally we searched for the positions D k R a    of the N virtual sources that minimized: Observe that in (5) and in (6) we were using that the virtual sources k a  were far from the boundary. The problem is that when the image is not supported far from the boundary this fitting function can create artifacts. To solve this artifacts arising from the virtual sources near the boundary we correct the fitting (7) changing S  by 1 11 That is, if we make a translation of the PSF to k a  , then m k is the mean of the translated PSF over all the pixels of the image (tile). In general we have 0  m k  1 / n and this number depends on each virtual source. When the source is inside the frame and far from the boundary m k = 1 / n, this amount decreases as the source approaches the boundary and m k  0 when the source is far outside. Now the function  2 is redefined changing in (7) S  by S . That is, the problem now is to find the N virtual sources D k R a    such that minimizes: Here we are using that all the random variables (x i ) are independent and have the same distribution. Then the mean over all the pixels it is equal to 1 1 Also we are using that the background is constant. As in SUPPOSe's method with B  0, we do not have a priori which is the total flux of the image (without background). We also have here the additional problem that we may have many virtual sources near the boundary. So we cannot use the relation of equation (5) in the first step of Algorithm 1. To replace that step, we develop a new sub algorithm to search for a relation between N and .
Sub Algorithm 1. We start with an initial  0 and T = S dev . At each step i, i. Calculate maxT and b k the point where it attains the maximum and ( ) iii. t (k) = || T ||. iv. The algorithm stops when t arrives to a minimum, and select N= Total iterations.
At the end With these new fitting function we can allow to implement SUPPOSe on images that are supported near the boundary or to make tiles without creating artifacts.

Results
In this section we will compare the results obtained after processing synthetic images with information near the boundary using SUPPOSe method and modified SUPPOSe method.
We also make a comparison of both methods when tiling real images acquired with objectives of different numerical apertures.

Reconstruction of synthetic images -SUPPOSe vs modified SUPPOSe
To test the method we start with the image synthesized in [12] consisting in two parallel straight segments 144 apart convolved with a Gaussian PSF, with noise and background added. The parameters used to synthesizethe image were selected to simulate the experimental conditions of a wide-field fluorescence microscope. The pixel size is 68 and the PSF half width (standard deviation) is 97.6. This is similar to the resolution obtained with emittedlight at a wavelength of 520 and a microscope objective of numerical aperture 1.3. In fig. 1b, the SUPPOSe reconstruction for N = 400 is presented. But, if we cut this image in two sub-images where on each image we use N = 200 we obtain the results of fig. 1c. We can see that the artifactson the boundary here are catastrophic. Instead when we apply the modified SUPPOSe method we obtain the results of fig. 1d. Many biological structure form a periodical arrangement. For example, fluorescence nanoscopy (STORM) recently revealed that actin, spectrin and accompanying proteins in axons form a membrane-associated periodic skeleton (MPS) arranged as periodic actin rings separated by  180 nm -190 nm, [49,50]. Nevertheless, as described in [50] this structure is not always present in all axons. This type of structures will cross the frame boundary and would present artifacts with the previous version for SUPPOSe. This motivated us to test our method with a synthetic family with a periodic structure. The purpose with this test is to prove that we can remove the boundary artifact and recover the period using the modified SUP-POSe. We built a family of synthetic images with 500 nm diameter axonal segments. Synthetic axons with presence of MPS structure are made up of point sources arranged in k parallel lines separated by a distance p, each line contains 50 point sources. This ground truth structure was convolved with a PSF Gaussian function to generate the synthetic image. The PSF Gaussian function model experimental conditions of a wide-field fluorescence microscope, the width of the PSF (standard deviation) is 95 nm and the pixel size is 34 nm (see Fig. 2a) and then was cropped to obtain an image with information close enough to the boundary (see Fig. 2b).
In this case the synthetic image was generated with k = 10 and p = 210 nm (see Fig. 2a) and then was cropped to obtain an image with information close enough to the boundary (see Figure Fig. 2b). In this example the real sources close to the boundary contribute with significant signal inside the frame. Observe that the total flux recorded by the image within the fitted frame is built by virtual sources well inside the frame, virtual sources inside the frame but close to the boundary (part of their flux falls outside the frame) and virtual sources outside the frame close to the boundary (part of their flux falls inside the frame). Fig. 1. (a)

b) Region of interest corresponding to the dashed white box in (a). (c) After applying SUPPOSe to the image in (b). We can see some artifacts near the boundary (dashed white box). (d) After applying the modified SUPPOSe to the image in (b). Scale bars 200 nm. Solutions convolved with a gaussian function of sigma 34 nm
After processing this part of the synthetic image with SUPPOSe we achieved the result shown in the Fig. 2c and observed some artifacts near the boundary of the image. If we apply the modified SUPPOSe to the same image, we obtain the result shown in the Fig. 2d.We show in this case that we can remove the boundary artifact present on the previous method and that we can recover the periodicity.
As we mentioned before, since the algorithm is parallelizable, for some large images it is useful to divide it into many tiles. Therefore, we want to test the accuracy of the modified SUPPOSe method when we increase the number of tiles. To this end, we again divide the synthetic image Fig. 2b into three tiles and apply the method of (10) . The results obtained using N  300 for each tile are shown in the Fig. 3. To evaluate the accuracy of the method, the positions of the virtual sources were projected in the direction of the synthetic lines on an axis directed along the axon long axis. Then we compute a histogram of this data set to see how they are grouped around the position of the ground truth lines. Finally we fit the histogram as a sum of Gaussian functions to show a multi-lobe distribution with averages for each lobe departing less than 1:7 nm from the ground truth, and standard deviations of 24 nm. We observe here how the artifacts disappear when tiling this periodic synthesized structure and that we also recover the periodicity.

Microscopy Images -SUPPOSe vs modified SUPPOSe
To test the performance of both methods on real images, we split a wide-field image of mitochondrial networks into two tiles and processed it with SUPPOSe and modified SUPPOSe. The image was acquired using a 10× objective with a low numerical aperture of 0.3 (Fig. 4a) giving a PSF half width (measured as standard deviation) of 420 nm of the objective for 599 nm light. The same region of the sample was captured using a 40× objective with a larger numerical aperture of 0.95 (Fig. 4b), in this case the PSF half width is 133 nm at the same wavelength. In this way, we obtained a ground truth structure to compare it with the SUPPOSe solutions. The solution obtained after processing the 10× images with the SUP-POSe is shown in Fig. 4c. The boundary artifact distorts the real structure, it can be clearly seen that in the area where both tiles meet, there is a break in the continuity of the mitochondrial network. Instead, when we process the same images with modified SUPPOSE these boundary artifacts disappear (Fig. 4d). Observe that the solution was convolved with a gaussian function of 133 nm standard deviation for its correct visualization, since it is the PSF width of the 40× objective. In this way we observe that we can resolve some structures such as the ring in the upper left area that cannot be resolved with the previous method, see Fig. 5. With this representation, we can conclude that the modified SUPPOSe has a resolution at least the same as the one of the objective of larger aperture. It is important to mention that the structures on the lower area of the frame are recovered with lower resolution since the signal to noise ratio is smaller in that zone. To show similar improvements in images taken with high NA objective we processed images of F-actin present in the endothelial cells of the bovine pulmonary artery (BPAEC) using NA = 1:3. Fig. 6 shows the results obtained from processing a specific region of the sample. We split the original image into two tiles and processed them separately with the two SUPPOSe methods. The results obtained are shown in the Figs. 6b and 6c. Artifacts generated by the original SUPPOSe method at the boundaries of the tiles distort information present in the whole region, as can be seen in the regions middle and left boxes in the Figs. 6b and 6c. Observe for example, that in Fig. 6b the region corresponding to the dash middle box may be confused with a structure of two straight segments. With the modified SUPPOSe the presence of a different type of structure in that region can be appreciated. Instead in the right box a structure that appears to be an boundary artifact is found that with the modified method it can be ascertained that it is the actual structure.

Conclusions
In this work we introduced a new fitting function for the method SUPPOSe presented in [12] to remove artifacts arising from the boundary. This correction allows parallelizing the process by tiling the image to improve the processing speed maintaining a good resolution in the intersection of the subfigures.
We first validated our method using synthetic images showing that in these cases the method solves the problem of the boundary artifacts. These artifacts are removed maintaining the same resolution that we had obtained when applying SUPPOSe to images supported inside the image frame. Finally we tested our method for real imag-es obtained from a standard fluorescence microscope. We could distinguish structures that were distorted by the previous method in F-actin filaments. We also remove artifacts and reconstruct structures of mitochondrial networks obtained by a low numerical aperture 10× objective. In this case we use the images obtained by a higher numerical aperture 40× objective as ground truth to compare showing thatthe method at least equaled the resolution of the high NA objective.