A joint Richardson—Lucy deconvolution algorithm for the reconstruction of multifocal structured illumination microscopy data

We demonstrate the reconstruction of images obtained by multifocal structured illumination microscopy, MSIM, using a joint Richardson–Lucy, jRL-MSIM, deconvolution algorithm, which is based on an underlying widefield image-formation model. The method is efficient in the suppression of out-of-focus light and greatly improves image contrast and resolution. Furthermore, it is particularly well suited for the processing of noise corrupted data. The principle is verified on simulated as well as experimental data and a comparison of the jRL-MSIM approach with the standard reconstruction procedure, which is based on image scanning microscopy, ISM, is made. Our algorithm is efficient and freely available in a user friendly software package.


Introduction
Optical super-resolution, the resolution of microscopic details below the diffraction limit, is revolutionising biomedical research. Since the first paper on stimulated emission depletion microscopy, STED, was published [1], a multitude of super-resolution techniques have arisen, all capable of 'breaking the diffraction limit' [2][3][4][5][6][7][8][9][10]. One of the most prominent methods in this field is structured illumination microscopy, SIM [7]. In its traditional implementation a sinusoidally modulated intensity pattern is used for sample illumination. Spatial frequencies in sample and excitation structures mix, and high frequency details contained in the sample become encoded in low frequency beat patterns, from which a high resolution image can be reconstructed. A particular modality of structured illumination microscopy is called multifocal SIM, MSIM [11]. In this variant multiple diffraction limited spots, rather than sinusoidal illumination patterns, are used to achieve super-resolution and the pixels on a widefield detector act as tiny pinholes, akin to a parallelised version of confocal microscopy. In fact, MSIM is the parallelised version of image scanning (confocal) microscopy, ISM [12], which uses the inherent super-resolving potential of confocal microscopy to double image resolution, when the detection pinhole is nearly closed [13]. For these reasons, the conventional algorithm for MSIM reconstruction [11] uses a confocal image formation model: where i x ( ) ⃗ is the recorded image, s x ( ) ⃗ is the fluorophore structure, e x ( ) ⃗ is the excitation light, and x PSF ( ) det ⃗ the detection point spread function, which limits resolution due to diffraction by the objective lens. The symbols · and ⊗ denote multiplication and convolution operations, respectively. In the case of ISM the excitation pattern is a single diffractionlimited spot such that = e x ( ) PSF ex ⃗ , where PSF ex is the illumination point spread function. In MSIM a regular pattern of such spots is used to parallelise the process. The conventional reconstruction algorithm in MSIM disassembles the recorded data from a parallelised to a sequential form and processes every single illumination spot individually to restore a super-resolved image. It must be noted, however, that equation (1) is only applicable in the case of an almost closed detection pinhole, while in general the widefield acquisition model has to be used. In this model the detected signal is the product of the sample structure and the illumination A joint Richardson-Lucy deconvolution algorithm for the reconstruction of multifocal structured illumination microscopy data pattern, which is convolved with the detection PSF. Using this model, we present here a joint Richardson-Lucy deconvolution algorithm [14], jRL-MSIM, which offers superior performance in reconstructing noise corrupted data compared to the standard ISM based method. Our MSIM reconstruction model treats all illumination points together as a widefield illumination structure. Hence, we demonstrate that there is a close conceptual analogy of MSIM with conventional SIM, which uses sinusoidal illumination patterns. We show that, as is the case for SIM, noise is efficiently suppressed also in our MSIM algorithm using joint Richardson-Lucy deconvolution. In [14], Ingaramo et al demonstrated the reconstruction of ISM images with jRL deconvolution and proposed that MSIM data could be reconstructed in the same manner by decomposition of the raw MSIM data to small areas containing individual spots. Each area can then be treated like a raw ISM image, and successive processing and combination of all spots, one by one, then generates the super-resolved image. In the approach we show here, we keep the high performance to deal with image noise but, because the algorithm uses a widefield model, it is inherently parallel and computational steps required are correspondingly reduced. In short we present an MSIM reconstruction procedure based on jRL deconvolution that is capable of processing the raw MSIM data without breaking it down into individual spots. We start by briefly reviewing the theory of realising super-resolution in confocal microscopy and ISM, as well as in SIM, and continue by reviewing the widefield models of SIM and MSIM and their similarities. Then we derive the concept of the widefield jRL reconstruction procedure, evaluate it on simulated as well as experimental data, and compare it to images reconstructed with previous MSIM algorithms.

Methods
For a well focused excitation beam in confocal microscopy the illumination spot becomes diffraction limited and we write = e x ( ) PSF ex ⃗ , where PSF ex is the excitation PSF. According to equation (1) it is possible to double resolution in a confocal microscope when the detection pinhole in the confocal plane is almost closed. Only then does the detection process become purely diffraction limited and, since excitation and emission light are passed through the same objective, PSF ex ∼ PSF det , i.e. the excitation and detection PSFs become comparable. Hence equation (1) becomes Equation (4)  their convolution will also be a finite function but with a supported frequency region whose cut-off k c conf is defined by the sum of k c ex and k c det , i.e. effectively doubled since k k c c ex det for a very small pinhole (see figures 2(a)-(c) for a comparison of imaging with an open and an almost closed pinhole). Hence the resolution is improved. As can be seen in figure 1, high frequencies near the superresolution cut-off k c conf are much less supported than the respective high frequencies in widefield imaging, k c wf , and are consequently much more affected by noise, especially Poisson noise. Use of such a small pinhole, however, leads to photon loss and hence reduces the signal-to-noise ratio, SNR, severely. For super-resolution microscopy this becomes prohibitive as seen on simulated data shown in figure 2, in particular figure panels d and e. In the ISM implementation the SNR is increased compared to a confocal set-up through the use of a two-dimensional array detector, e.g. a CCD, instead of a photodiode situated behind a physical pinhole. Hence, each pixel on the camera now acts as an individual, almost closed, pinhole. The pixels that lie off the optical axis with a shift o ⃗ detect the same information as the central, confocal, pixel, but with the signal shifted by o / 2 ⃗ as seen in figures 2(f)-(h) (depending on the Stokes shift between excitation and emission wavelengths, the exact off-set may be smaller [16]). Since many adjacent pixels gather signal simultaneously, there is no inherent loss in photons, and hence conf of the OTF conf is twice as large as that of OTF wf , hence doubled resolution can be achieved. Frequencies near the respective cut-offs, however, are much less supported in the confocal case than in the widefield case and hence more affected by noise. SNR, with this technique. The image formation on every individual pixel is thus described by [12] ex det (5) Reassigning the information from off-axis pixels to positions half-way towards the on-axis pixel, i.e. shifting pixel signals by o / 2 ⃗ , is then equivalent to overlapping shifted but otherwise identical image structures. This increases the SNR dramatically. A detailed analysis of this approach is presented in [16]. Furthermore, since every single pixel acts as an almost closed pinhole one retains the resolution-doubling potential of confocal microscopy described in [13] and in the previous section, however without inherent loss of signal. Scanning the confocal spot across the whole field of view and repeating this reassignment procedure at all beam positions creates the final super-resolved image. MSIM and similar techniques [16][17][18] parallelise the ISM principle by sweeping hundreds of individual focal spots in parallel to achieve higher imaging speeds, however the raw images are then processed for each confocal spot according to the ISM principle above.
Viewing MSIM as a patterned widefield imaging technique instead opens the door to a different reconstruction approach. In conventional SIM a sinusoidal illumination pattern, ( ) k x cos ex ⃗ ⃗ , is used to mix spatial frequencies of the excitation pattern with the underlying structure [7]. The widefield model becomes: ex det (6) in real space, and: in Fourier space. Constant factors were omitted in equation (7), and k ex ⃗ denotes the wave-vector of the sinusoidal illumination pattern. splitting the acquired raw data into their individual frequency bands and reassigning them to their proper locations, ± k ex ⃗ , in frequency space. It is the shifting operation of delta pulses in Fourier space that enables the resolution enhancement. The Fourier representation of multifocal illumination patterns as used in MSIM can also be disassembled into a multitude of delta pulses by treating the illumination pattern as a diffraction-limited delta comb, Ш T , with period T: The Fourier transform of Ш T is Ш T 1/ , again a delta comb but with periodicity of 1/T, which describes the delta pulse spacing in frequency space. To split the individual frequency bands, a number of individual images need to be acquired with spot patterns appropriately shifted between acquisitions. The individual split frequency bands can then be repositioned and again the resolution is increased twofold. However, because a much greater number of frequency components have to be extracted and repositioned, this operation is far more complex than in the standard SIM case.
A sophisticated way of processing widefield SIM data was introduced by Ingaramo et al [14] in the form of joint Richardson-Lucy deconvolution, jRL. The algorithm performs the reconstruction without the need of band-separation and repositioning and incorporates an image deconvolution process. Furthermore, as jRL is an improved version of the RL deconvolution algorithm [15], it is particularly capable of handling images corrupted by Poisson noise. The jRL algorithm requires a mathematical model of the image formation process and precise knowledge of the illumination structure to iteratively combine a multitude of single raw data files into the final super-resolved image. It maximises the likelihood of an estimated super-resolved image ŝ to be the source of the measured raw images i when imaged in the presence of Poisson noise: As shown above MSIM data generated with spot patterns can be represented in a similar manner as conventional SIM, where sinusoidal patterns are used for illumination. Hence, we use the MSIM widefield model with delta comb patterns as input for Ex m . To reject out-of-focus light and allow optical sectioning we additionally pinhole the raw data by multiplication with Gaussian shaped digital pinholes before starting the actual reconstruction process. Each step of the iteration provides an update of the estimate s x ( ) ⃗ that is a closer approximation of the original structure s x ( ) ⃗ :  [14]. Our algorithm, which we subsequently refer to as jRL-MSIM, then combines in a natural fashion a widefield model of MSIM with a deconvolution procedure to reconstruct super-resolved images.

Results and discussion
As a proof of principle we created artificial MSIM data using the presented widefield model and processed it with the jRL-MSIM reconstruction algorithm which we implemented in MATLAB. The simulation assumed an illumination wavelength of 488 nm and a pixel size of 40 nm.
The resulting image, displayed in figure 3(a) right panel, shows a clear resolution improvement in comparison to the respective widefield image, figure 3(a) left panel, and demonstrates the validity of the model for processing. In this case the parameters of the illumination pattern were known exactly and noise free raw images were used. In figure 3(b) an enlarged region of the jRL-MSIM result is furthermore compared to widefield and ISM based [11] reconstructed images; the original image is also shown. For a valid comparison of the jRL-MSIM reconstruction to both the widefield and ISM reconstructed images, the latter images shown in figure 3(b) were deconvolved with a conventional Richardson-Lucy deconvolution algorithm. The jRL-MSIM image visually appears to be superior to the ISM based reconstruction and comes closest to the original image, but clearly both ISM and jRL-MSIM algorithms produce super-resolved images that are significantly sharper than the widefield case. The slightly inferior ISM performance may be reflective of accumulated rounding errors in the pointwise processing procedure. In practise, however, neither the illumination To capture more signal and improve the SNR it is possible to use multiple off-axis pinholes and take several images. As all of these images are captured using 'closed' pinholes they will have the same high resolution but are off-set from one another (f). Summation of a multitude of off-set 'closed'-pinhole images results in an almost noise-free but non-super-resolved image (g). Repositioning the raw images before their summation, however, maintains the high resolution while simultaneously increases the SNR (h). For an excellent and more detailed overview of the fundamentals of MSIM we refer the reader to work by York et al [11,18].
parameters are known precisely nor is the recorded raw data free of noise. Consequently, to test the jRL-MSIM algorithm on such cases, we added Poisson noise to the raw data with an expectation value of 128 photons in the brightest pixel and let the software estimate the illumination parameters directly from the raw data using the algorithm presented by York et al [11]. Again, we compared the reconstructed super-resolved images created by the two reconstruction methods. The results are shown in figure 3(c) and demonstrate clearly the superiority of jRL-MSIM in this case.
In some cases, the illumination parameters cannot be estimated with sufficient precision and this leads to artefacts in ISM based reconstruction algorithms. This is shown in figure 4, which compares the two methods for identical illumination parameters. The ISM reconstructed image suffers from overlaid pattern noise, which is not present on the jRL-MSIM reconstruction. Instead, the jRL-MSIM image appears crisper and free of artefacts.
For a comparison of computational efficiency the actual illumination pattern has to be taken into consideration. For 143 raw images, each featuring 256 × 256 pixels with ca 100 illumination spots, the reconstruction and successive deconvolution using the ISM algorithm took 25 s, compared to 5.6 s for a single iteration of jRL-MSIM. We tested jRL-MSIM on real images featuring different SNRs and found that optimal results were obtained for between 5-30 iterations in all cases. The computational efficiency of jRL-MSIM becomes superior to that of ISM when the spot density in the illumination pattern is increased. Here the speed of jRL-MSIM increases since its processing time scales only with the number of images while that of ISM-MSIM scales with the total number of processed spots. In addition, the smaller number of raw images required reduces the acquisition time of the data sets. Note, that the spot density must remain below a value where signal crosstalk between corresponding illumination spots occurs.
To validate our reconstruction approach in practise, we applied the jRL algorithm to experimental image data obtained of fluorescent beads and biological samples. The MSIM raw images were obtained in a home built MSIM microscope similar to the one described in [11] using a spatial light modulator for pattern generation and laser excitation of 488 nm wavelength. Figure 5(a) shows the algorithm's performance on 170 nm in diameter fluorescent beads stained with Alexa Fluor 488 dye. Clearly the resolution is enhanced by jRL-MSIM. The first tile in figure 5(b) shows a widefield image created from publicly available MSIM data of labelled microtubules provided by Andrew York (code.google.com/p/msim/). The second tile shows reconstruction results using the point-wise ISM reconstruction procedure presented in [11] and the third the results of jRL-MSIM using an underlying widefield model. Note that again both algorithms perform well in enhancing the resolution.
The algorithm is capable of analysing threedimensional data as well. Although jRL-MSIM does not explicitly improve the resolution along the imaging axis, digital pinholing leads to a significant contrast enhancement along z, which is a property well known in confocal microscopy. The rejection of out-of-focus light reduces blur and therefore causes an apparent increase in resolution. This is shown in figure 6. The figure shows 3D renderings of two human red blood cells, imaged in MSIM mode. The cells were stained following standard protocols. The left rendering shows the sum of all raw images without processing, akin to what would be obtained in taking an image stack in widefield mode, and the right panel shows the jRL processed image stack. Clearly the overall appearance is much improved in resolution for the latter. Furthermore, from the cross sections through the rendered volumes a resolution improvement along the optical imaging axis is clearly observable: the overlap and relative positioning of the two red blood cells is clearly discerned in the jRL-MSIM processed data due to the much improved contrast obtained.

Conclusion
Both processing techniques, the confocal ISM based MSIM reconstruction algorithm as well as the widefield  (a) 170 nm beads stained with Alexa Fluor 488 were imaged using multifocal excitation-patterns and processed using jRL-MSIM. The jRL image achieves twofold resolution compared to the respective, deconvolved, widefield image. (b) Microtubules stained with with Alexa Fluor 488 dye were reconstructed using the standard multifocal ISM reconstruction procedure [11] and jRL-MSIM. Both techniques achieve comparable super-resolution and successfully reject out-of-focus light. jRL-MSIM approach, succeed in out-of-focus light reduction and they achieve similar performance to enhance resolution when reconstructing noiseless images. In the presence of noise, however, the estimation of illumination parameters is more challenging and the jRL-MSIM algorithm proves to be superior. Additionally, depending on the number of iterations, the jRL-MSIM algorithm can be faster than ISM-MSIM. Finally we demonstrate the close conceptual relationship of the MSIM spot scanning technique to both the widefield sinusoidal SIM and the confocal ISM techniques and show an alternative way to achieve super-resolution and optical sectioning using multifocal illumination patterns. The reconstruction software which implements both ISM-MSIM and jRL-MSIM algorithms in an easy to use software package, including manual and example files, is freely available at laser.ceb.cam.ac.uk. Figure 6. 3D renderings of two human red blood cells, imaged in MSIM mode. The left rendering shows the sum of all raw images to produce a widefield image; the right rendering was generated using the jRL-MSIM algorithm. Especially in the cross sections through the volume the successful rejection of out-of-focus light and higher contrast through digital pinholing and the additional apparent increase in resolution in the z direction can be seen.