Pupil mask diversity for image correction in microscopy

: Three-dimensional microscopy suﬀers from sample-induced aberrations that reduce the resolution and lead to misinterpretations of the object distribution. In this paper, the resolution of a three-dimensional ﬂuorescent microscope is signiﬁcantly improved by introducing an amplitude diversity in the form of a binary amplitude mask positioned in several diﬀerent orientations within the pupil, followed by computer processing of the diversity images. The method has proved to be fast, easy to implement, and cost-eﬀective in high-resolution imaging of casper ﬂi:GFP zebraﬁsh.


Introduction
Imaging fluorophores inside large three-dimensional samples with sub-cellular resolution is of a great interest in life sciences, however the inherent optical inhomogeneity of biological tissues distorts the image. Aberrations, arising from spatially variant phase delays throughout the sample, reduce the image contrast and scramble both the excitation and the fluorescent light, resulting in images that misrepresent the real fluorophore distribution.
These misrepresentations, in the fluorescent imaging arm, may be corrected only by first elucidating the error and then either by physically compensating for the phase delay using adaptive optics (AO) [1] or by computationally reconstructing the distribution, known as deconvolution [2].
Whilst light-sheet fluorescence microscopy (LSFM) [3,4] is a technique that has been specifically developed for imaging large three-dimensional samples, the aberrations, scattering and absorption only generally allow for the use of lower numerical apertures (NA) that dampen the effect of aberrations at the cost of resolution. Despite this, sample-induced fallaciousness is present in unprocessed LSFM image stacks and researchers tend to acquire a number of datasets from different "views" and computationally merge them together [5]. To further improve the image quality these datasets may also be deconvolved with the measured system point-spread functions (PSFs) [6].
Adaptive optics (AO) has been employed more recently to correct for these aberrations in both detection [7] and excitation [8], where it proceeds by a process of wavefront sensing and correction. In microscopy specifically, there has been two methodologies proposed for the sensing of aberrations. The first directs light to a sensor, such as a Shack-Hartmann sensor [9]. This approach requires the use of a "guide star" [10] and is better suited to scanning techniques rather than wide-field. The second method is an indirect measurement via the optimisation of the image quality and is preferred in wide-field imaging because using a wavefront sensor on an extended source is problematic.
Image optimisation requires the quantification of image quality in the form a metric. The choice of such a metric in wide-field microscopy is problematic and often sample and technique dependent. The optimisation of the adaptive optical element's (AOE) control inputs, therefore, does not always lead to the image with optimal fidelity to the fluorphore distribution. The reasons for this is that in these microscopes the light source is never truly confined to a single-plane and the aberration is spatially variant.
These two points imply that the light received at the camera is a superposition of light from different regions within three-dimensional sample that would require different control inputs to correct them. As a result, it would take at least N + 1 frames [11] where N is the number of degrees of freedom of the AOE, to compensate the image in a single region of the sample with a specific aberration. In the lingo of AO these regions are called isoplanatic patches.
Given the effect of photo-bleaching and photo-toxicity, it is often better to perform deconvolution on such wide-field images to provide a better estimation of the object when the point spread function (PSF) is known even approximately. In the presence of unknown aberrations, however, the PSF is unknown and deconvolution with the incorrect PSF will result in a poorly estimated object and unintended artefacts. A blind single-frame deconvolution [12] is an ill-conditioned mathematical problem with multiple solutions and therefore, cannot be done reliably in the presence of unknown aberrations without substantial a priori information. As a result, deconvolving an LSFM image stacks with a measured system PSF will not yield the true fluorophore distribution because of the presence of aberrations.
In order to deconvolve in complex three-dimensional samples correctly, with aberrations, it is necessary to either: (1) correctly identify the point-spread function by using diversity -this produces a multi-frame blind deconvolution problem [13] that is better conditioned and solutions can be more reliably found via a plethora of bilinear optimization algorithms [14]; (2) or correct for the aberrations using adaptive optics such that the PSF is known to be a diffraction-limited Airy pattern before deconvolving.
In this article, we develop an extremely simple and efficient implementation of the first method. By combining physical perturbation of the optical system, done by a rotating amplitude mask, alongside the computational processing of the resultant images, an image stack is corrected for sample-induced aberrations. To the authors' knowledge, whilst phase diversity has been used for astronomy [15], horizontal imaging [16], ophthalmic imaging [17], to calibrate deformable mirrors [18], and with amplitude diversity in the wavefront sensing [19]; there has been no successful attempt to apply this general sort of diversity in microscopy. The reasons for this may be that it has been difficult to identify the point-spread functions and obtain a well-conditioned deconvolution without noise-amplification at low signal-to-noise ratios (SNR).
An outline of the article: first an overview of the theory is given, then the description of the experimental methodology, the results from multiple casper fli:GFP zebrafish imaging experiments and to end a short discussion of the approach.

Methodology for aberration correction
An understanding of the technique may be derived from elementary imaging theory. In an adaptive optics system with control inputs x, the relationship between the formed image i(x) and the pupil function P(x) = A(x) exp ( jφ(x)), where A(x) is the aperture's amplitude function and φ(x) its phase function, may be modelled by the following convolution operation: where w is the noise and o is the unknown fluorphore distribution. Normally, the goal in AO is to find x * such that one has i(x * ) = | F{A(x * )}| 2 * o -the diffraction-limited image. From an imaging point of view, it may be noted that x * itself is not important, but instead the fluorophore distribution o. It is possible to find an estimate of o without finding x * by solving a mathematical optimisation problem. For control inputs x, the following problem may be posed: where only i(x) is known and h(x) = | F{P(x)}| 2 . This is a many-to-one problem since there are many h(x) and o that may produce the same i(x). Its inverse is an ill-posed problem; however, if enough is known a priori one can solve this problem and crucially in an AO system one can acquire images with different x, adding diversity and therefore, reducing the number of solutions that satisfy the equation.
To understand why this is necessary, one must consider that in the presence of aberrations, the OTF H(x) = F -1 {h(x)} has points which are heavily modulated, implying that the system no longer transmits a set of spatial frequencies {ν ϕ }. In this case, I(x) = F -1 {i(x)} contains no information about the object O = F -1 {o} at these frequencies and no computational process can hereby restore this information. Furthermore, all image registration processes have additive noise both due to the sensor and the quantum nature of light, which is clearly prevalent in low light applications such as fluorescence microscopy. The image spectrum I(x) may then be regarded as the sum of the true image spectrum I r (x) and the noise spectrum W = F -1 {w}: in these spectral regions, destructive noise-amplification occurs as can be seen by: Clearly, very large or infinite frequency components are non-physical and come from the inverse relationship -by adding the perturbations to the aperture one may move these zeros and extra information about the object may be included in the estimate by using multiple diverse observations.
In the notation the x shall be dropped and a subscript m used to denote a frame acquired under a different x value. The process begins by acquiring the first frame I 1 , which also becomes the first estimate for the objectÔ (1) . For all subsequent frames acquired within the iterative loop, the k th step for the m th image, where 1 ≤ m ≤ k, the Fourier relationship is: The estimate of the OTFsĤ (k) m are updated using the previous estimate for the object spectrum where P H is a projection operator onto the space of real-valued non-negative PSFs of a particular spatial extent. With this new estimate of the OTFs an estimate for the object spectrum may be obtained by back-projection using a window of the last n measurements: where 0 < α < 1 is a feedback parameter. The PSF and thus the image is perturbed in the next step and this allows the algorithm to begin with a realistic starting point for O and H. The choice of Eq. (7) and Eq. (6) are due to noise propagation and discussion of the particular form of operators P H and P O are beyond the scope of this article, but may be found in Wilding et al. [20]. The form of the algorithm has been altered from Wilding et al. [20] to run recursively on new data, hereby, it has been called the recursive Tangential Iterative Projections (rTIP) algorithm in this use. For this reason a more systematic representation has been given in Table 1. Here the steps of the algorithm along with concise descriptions have been recorded. The feedback parameter's function in Eq. 7 is to bias the solution towards the new estimate (α > 0.5) or the previous one (α < 0.5). Depending on the window size, one may be more certain of the previous solution, based on many frames, than the new one based on only a couple of measurements. It need not be static, as shown in Step 5 in the Table 1, where α is a measure of the average relative difference between the new initial estimateÕ (k) and the lastÔ (k−1) . This term provides a penalty for new estimates that differ from the previous ones and slows down the process of aberration correction, but reduces noise amplification. Table 1. Steps of the rTIP algorithm with their descriptions in schematic form. After the second acquisition the steps 1 to 5 are repeated with every subsequent acquisition to update the object estimateÔ (k) .
Step/Projection Process Description Step Step Multi-frame linear deconvolution Step Inverse Fourier Transform Step Step

Experimental design
A light-sheet fluorescence microscope has been realized in the configuration shown in Fig. 1, the excitation light is provided by a 488nm laser (100mW Sapphire LP, Coherent Inc., U.S.) and polarization optics (WPH10M-488 and GT5-A, Thorlabs, U.S.) to produce linearly polarized light for the spatial light modulator (SLM) (512x512, Meadowlark Optics, U.S.). The beam of the laser is split to multiple setups and therefore, is run giving 9mW into the LSFM microscope path, where with the efficiency of the SLM due to the grating width gives around 1mW at the back aperture of the objective. The SLM is conjugated to the back aperture of an NA= 0.3 objective lens (UMPLFLN 10x Olympus, Japan) via a beam expander (AC508-200-A-ML and AC508-200-A-ML, Thorlabs, U.S.). The light-sheet is formed by using a cylindrical lens pattern applied to the SLM that reduces the focal power of the objective in one direction. Orthogonal detection of the fluorescence light is done with a NA= 0.5 imaging objective (UMPLFLN 20x Olympus, Japan), via the rotating actuator (PRM1/MZ8/TDC001, Thorlabs, U.S.) close to the back aperture with a 3D printed acrylonitrile butadiene styrene (ABS) double-wedge (Fig. 2(d To design the pupil mask simulations of many different types of pupil masks where tried from random binary masks to spirals and wedges. The wedge design was chosen for its following properties: it was easy to manufacture and mount; in a 360 • rotation of the mask all parts of the pupil would be transmitted; the diversity added by the mask empirically seemed to satisfy the condition of recovering lost frequencies m H m (ν ϕ ) > 0; and it was easy to change the diversity added by adjustment of the wedge angle. In this way, a set of wedges with different angles were produced and tested experimentally. The wedge with the best performance for a particular set of imaging conditions could then be chosen.
A wedge angle of 45 • was found to be the best pupil mask, balancing between diversity generated and light loss. The sample used was a para-formaldehyde fixed casper fli:GFP zebrafish (danio rerio) mounted in 0.5% agarose in a 1 × 1mm square capillary tube (Vitrotubes TM , VitroCom, U.S.) aligned with the LSFM axes. Three-dimensional acquisition is done by moving the sample through the light-sheet with a servo-motor stage (USB Stage, Picard Industries, U.S.) followed by post-processing registration to adjust for any sample drift.

Imaging results
A N = 8 acquisition proceeds by acquiring a single z-stack, rotating the wedge 45 • and then acquiring the next stack until the wedge has been rotated through the entire 360 • range. With a perfect double-wedge there would be degeneracy in the dataset, however, this is not the case with the design used due to imperfections in manufacturing and lack of perfect centering on the pupil around ±1mm. The resulting images at each z position have 8 frames with different diversities. These may be processed sequentially after three-dimensional registration using the methodology described in Section 2 and a resultant series of estimates based on increasing the number of frames may be yielded. Approximations of the aperture shapes used for each of the images is shown by Fig. 2(c).
By taking one of the planes (at z = 75µm depth) the restoration process can be demonstrated in Fig. 2(a), where the frames of the dataset, i.e. normal frames from the microscope are compared with the successive fluorophore estimates. In this figure a 256 × 128 px region-of-interest (ROI) of the original 2048 × 2048 px images is shown. It can be observed that the reconstruction sharpness and contrast improves with increasing number of frames, but not with the raw images.
A 360 × 198 px ROI at z = 135µm depth from the dataset without the mask is also shown in Fig. 2(b), where it can be seen that the introduction of the mask does, as expected, lead to a loss in image quality through loss in signal in the first masked frame, however, by continuing the acquisition process a better result is yielded. This improvement can be more clearly seen by taking the line profile (LP) shown in Fig. 2(e). The LP is plotted for no mask, frame 1 and frame 8 of the iterative processing. This plot highlights the recovery of high spatial frequencies and image contrast between the first and the last frame, and the appearance of features lost by the optical aberrations in the original LSFM image shown by the red arrow.
The restoration process may then be repeated either plane-by-plane, in smaller ROIs to combat anisoplanatic aberrations, or as a full three-dimensional deconvolution. All of these method will yield a final single three-dimensional dataset. Since the aberration is spatially variant in LSFM, the plane-by-plane method is preferred, allowing the aberration to change between each plane along the direction orthogonal to the plane. The PSF is expected to vary with depth z, but also with position across the light-sheet x, especially axially; therefore, x y processing and yz processing are tried.
The maximum intensity projection (MIP) of the unprocessed LSFM data and the processed set is shown in Fig. 3(a) showing the increased contrast of the processed datasets. Fig. 3(b) and (c) give an insight into the effect of the processing in the z-direction or the x-direction. It can be seen that the xy plane-by-plane processing does not improve the z-sectioning ability as expected, however, when processing the yz slices the axial resolution can be seen to be improved. This effect is most clearly seen in the line profile found in Fig. 3(c). Through these experiment it is shown that by adding a simple AOE it is possible to improve the images from a LSFM. This occurs because the active modulation of the OTF removes the zeros in the summed image spectrum ensuring that the object information may be extracted algorithmically. Unwanted noise amplification does not occur due to both the combination of the physical and computational process used (Section 2).

Discussion
The benefits of this new approach may be summarised at this point: firstly, it is not necessary to determine what are the control settings that optimise the image, this means that the amount of expertise required to use adaptive optics well is reduced dramatically; secondly, the method requires less frames thus incurring less photo-bleaching and toxicity than would be required by model-free or model-based optimisation of the images with most devices used for microscopy. Thirdly, it is possible to go beyond the limitations of the AOE, therefore, cheaper and less technical devices may be used for the same effect. Finally, the technique finds the spatially variant PSFs from the data and therefore, takes into account sample-induced aberrations instead of just the system ones, reducing the probability of yielding a solution that misrepresents the fluorophore distribution.
This is not to state that this technique is perfect, since one of the benefits of physically optimizing the PSF is that it improves the signal-to-noise ratio (SNR) of the images. Aberrations through the loss of SNR alone can cause a dramatic reduction in image quality especially in non-linear imaging modalities. It is therefore, a valid criticism that reducing light collection efficiency is to be avoided. The reason for making the trade-off here is that diversity can restore the lost spatial frequencies and therefore, aberrations can be corrected.
Furthermore, the speed of the algorithm is linearly proportional to the number of pixels used. For CPU it can process at 2.5µs per pixel (Intel i7 8-core, 32 GB RAM), meaning for the datasets in this paper of (2048 × 2048 × 100) px × 8 z-stacks, the total processing time is around 140 minutes. This is not an unreasonable for blind deconvolution on LSFM data, but slow for practical imaging. The problem has parallelism and the speed can be increased by an order of magnitude by processing on a GPU. For a single 2D plane of 8 images on the GPU (nVidia GTX980), it takes 6.9 seconds to have a corrected reconstruction, which is faster than performing image optimisation on the LSFM data in the author's experience.
The limitations to the correction depend on the signal-to-noise ratio and are difficult to quantify experimentally. Some numerical work has been on this question in the authors' previous work in Ref. [20] for the interested. One can imagine the aberration as providing another source of noise in the imaging process. As a result, the aberration size, its order, the brightness of the fluorophores and quality of the detection system all combine to limit the total dynamic range of the correction. As an intuitive example, a bright sample would be able to be corrected at a higher aberration level than a dimmer sample if it had the same SNR.
It should be noted that this type of pupil mask is not the only possible option. A transmissive element, such as a pyramid, could also be used in the case where photons are too precious to lose by collecting all the light paths. Or if a phase element is used such as a deformable mirror, the technique may also be used in parallel with optimisation-based approaches to obtain the corrected image faster. To do this, however, it is preferred to use an optimisation algorithm that does not rely on a coordinate search since for best results one requires the root mean square phase diversity (r.m.s.) to be less than the r.m.s. phase aberration, such an algorithm may be found in Verstraete et al. [21]. In this hybrid-case, the bonus is that the physical Strehl ratio and thus the SNR ratio improves, which may be useful in low-light applications and on light-sensitive samples.
To conclude, in this article a new method for applying adaptive optics in the detection path of a wide-field microscope has been presented for the purpose of obtaining a more accurate representation of the fluorophore distribution. It relies on the principle of changing the aperture between stack or image acquisitions, hereby changing the PSF of an optical system so that it is possible to computationally extract the object distribution through bilinear optimisation. The technique provides a faster, more cost effective solution to the problem of aberration correction imaging in wide-field microscopy when compared with existing methods, however, it does suffers from some of the common drawbacks of computationally-based adaptive optics such as there is no physical increase in the signal-to-noise ratio and the scaling of the problem with image size.