Cell imaging beyond the diffraction limit using sparse deconvolution spatial light interference microscopy.

We present an imaging method, dSLIM, that combines a novel deconvolution algorithm with spatial light interference microscopy (SLIM), to achieve 2.3x resolution enhancement with respect to the diffraction limit. By exploiting the sparsity of the phase images, which is prominent in many biological imaging applications, and modeling of the image formation via complex fields, the very fine structures can be recovered which were blurred by the optics. With experiments on SLIM images, we demonstrate that significant improvements in spatial resolution can be obtained by the proposed approach. Moreover, the resolution improvement leads to higher accuracy in monitoring dynamic activity over time. Experiments with primary brain cells, i.e. neurons and glial cells, reveal new subdiffraction structures and motions. This new information can be used for studying vesicle transport in neurons, which may shed light on dynamic cell functioning. Finally, the method is flexible to incorporate a wide range of image models for different applications and can be utilized for all imaging modalities acquiring complex field images.

1. Introduction

Classical light microscopy techniques cannot be used directly in imaging most biological structures, as they do not significantly absorb or scatter light [1]. Interference-based methods such as phase contrast [2] and differential interference contrast microscopy [3] allows imaging these transparent structures without the need for staining or tagging. Recently, more advanced methods introduced the ability to measure quantitative information on the specimen by precisely quantifying optical phase shifts induced by the structure and motion of the specimen [4]. Spatial light interference microscopy (SLIM) [5], is a new and powerful quantitative imaging technique which allows high phase sensitivity imaging of nanoscale structures. SLIM has the important advantages of utilizing illumination with short-coherence length, and the ease of implementation via add-on modules on existing phase-contrast microscopes.

Although interference-based microscopy has tremendous advantages, it is still affected by the optical degradation and noise introduced by the instrument [6]. These degradations can be removed to a certain extent by employing post-processing methods. Deconvolution is a common postprocessing method to invert the optical transfer function of the instrument. Although it is widely used in intensity-based microscopy [7-11], not much work has been reported on deconvolution in microscopy systems collecting quantitative information through complex fields. The work in [12] investigated the use of complex field deconvolution through inverse filtering in digital holographic microscopy [13], and have reported that the noise amplification, commonly encountered with inverse filtering in intensity-imaging, is not as significant in the case of complex field microscopy. A nonlinear deconvolution method has been developed in [14] for SLIM that estimates the unknown magnitude and phase fields via a combination of variable projection and quadratic regularization on the phase component.

In this paper, we present a novel method, dubbed dSLIM, for complex field deconvolution using an image model suitable for characterizing the fine scale structures. Based on the prominent features of phase images, we model the underlying image using the sparsity properties of the transform coefficients. This model is especially useful in capturing fine-scale structures, and it successfully reveals the details in the phase components lost due to the instruments optical transfer function. In addition, due to the very low noise floor provided by SLIM (0.3 nm spatially and 0.03 nm temporally [5]), accurate experimental estimates of the point spread function can be obtained, and deconvolution artifacts are significantly reduced. We demonstrate that resolution increases by a factor of 2.3 can be achieved with dSLIM. Thus, in dSLIM, images with a final resolution of 238nm can be rendered with only a 0.65NA objective. Moreover, the presented methodology is very flexible in incorporating image properties in a wide range of applications (such as materials imaging) and can also be utilized for other interference microscopy techniques.

The rest of this paper is organized as follows. We provide an overview of SLIM in Section 2. The image degradation model and the general framework for complex deconvolution is presented in Section 3. The proposed deconvolution algorithm dSLIM is developed in Section 4. We present experiments with SLIM images in Section 5 and conclude in Section 6.

We use the following notation throughout the paper: Bold letters h and H denote vectors and matrices, respectively, with transposes hT and HT. The spatial coordinates within a image are denoted by (x,y), operator * denotes convolution, and i is equal to . Finally, {·} is used to denote a set created with its argument.

2. Overview of Spatial Light Interference Microscopy (SLIM)

In SLIM [5], a spatially coherent light source, U(x,y) = |U(x,y)|exp[iΦ(x,y)], is used for illumination, which is decomposed into scattered and unscattered fields after passing the specimen. Let us denote the unscattered light as U0 and the scattered light as U1(x, y). A liquid crystal phase modulator (LCPM) is used to introduce phase modulations to the unscattered field, such that where ϕ0 is the intentionally added phase delay, and ϕ1(x, y) is the phase difference between the scattered and unscattered fields caused by the specimen. The unscattered field contains the uniform background of the image field, whereas the scattered light provides information on the structure of the specimen. The recorded intensity is expressed as In traditional phase-contrast microscopy [2], ϕ0 is fixed at and a single image is acquired. This can only provide qualitative information (i.e., ϕ1(x, y) cannot be uniquely retrieved). In contrast, SLIM uses multiple phase delays 0, , π, and , such that ϕ1 can be uniquely determined. Specifically, ϕ1 can be extracted from the four recordings using Moreover, the phase associated with the complex field can be calculated by where we define by the ratio of magnitudes of the scattered and unscattered light.

3. Image formation and deconvolution model

As demonstrated in the previous section, both the magnitude and phase of the complex image function U(x,y) = |U(x,y)|exp[iΦ(x,y)] can be uniquely determined using SLIM (using Eqs. (3) and (5)). However, as in all imaging systems, only a degraded version of this field can be observed in practice. Modeling the imaging process as a linear, spatially invariant degradation system, the measured image can be expressed as the convolution of the original complex field with the instrument point spread function (PSF) as where h(x,y) is the PSF of the system, and n(x,y) is the additive signal independent noise. In general, both the magnitude and phase of the complex image function is degraded via the optical transfer function. As in traditional deconvolution [10], these fields can be estimated using a regularized inverse formulation where σ2 is the noise variance, and the functional R(·) is used to regularize and impose constraints on the estimates of the magnitude and phase.

The estimation of both the magnitude and phase in Eq. (7) is a nonlinear optimization problem, which is highly ill-posed and hard to solve in practice. However, in practice, the degradation in the magnitude is very small compared to the degradation in the phase [12]. In addition, the phase image contains most of the information of the specimen relevant and useful from an application point of view. Therefore, it is convenient to assume that the magnitude of the field is constant, and the imaging system introduces negligible distortion in the magnitude, such that |Ũ (x,y)| ≈ |U(x,y)| = const. This assumption makes the problem in Eq. (7) linear, and it is also very useful in avoiding instabilities due to nonlinearity. With this approximation, the problem Eq. (7) becomes For mathematical convenience and clarity, let us denote by g(x,y) the observed field exp[iΦ̃(x,y)], and by f (x,y) the unknown field exp[iΦ(x,y)]. Due to the linearity of the degradation, the problem Eq. (8) can be expressed equivalently in matrix vector form as where g and f are images g(x,y) and f (x,y) in vector forms, respectively, and H is the convolution matrix corresponding to the PSF h(x,y).

When no regularization is used in Eq. (9), the closed-form solution of f̂ can be found as (HTH)−1 HTg (equivalent to the inverse filter). However, this approach generally leads to noise amplification and ringing artifacts due to heavy suppression of high spatial frequencies. The role of regularization is to impose desired characteristics on the image estimates to avoid this noise amplification and to increase the resolution. The parameter β is used to control the trade-off between the data-fidelity and the smoothness of the estimates.

4. Complex field deconvolution using sparsity

4.1. Image model

It is well-known that phase-contrast imaging is highly sensitive to object boundaries but relatively insensitive to the flat background areas. Due to this characteristic, phase images exhibit high contrast around edges and local spatial variations within the specimen, which allows the capture of accurate shape and edge information. Moreover, in many important applications where interference microscopy is utilized, e.g., live cell imaging, the specimen contains very fine structures and small-scale movements that need to be precisely localized.

Based on these observations, our goal is to construct a model of phase images that accurately represents this fine structure with sharp boundaries. We base our modeling on the sparsity principle, that is, our main assumption is that the phase images can be very accurately represented in some transform domain with sparse coefficients [17, 18]. This transform sparsity can be achieved by appropriately selecting the transforms that capture the characteristics of spatial variations within the image.

In this work, we consider a set of L linear transforms D of the complex image f with k = 1,...,L. These transforms are chosen to be high-pass filters, such that their application provide complex images with a large number of coefficients with small values with only a few coefficients containing the most of the signal energy. The selection of the linear transforms that most accurately capture the image characteristics is crucial in the final image quality. We employ a collection of difference operators to capture signal variation at varying scales. The directional contourlets [19], depicted in Fig. 1, are used to capture the overall spatial variation, as they contain both vertical/horizontal and diagonal directions. In addition, for smaller scale features, we include first order directional difference operators and 45° and −45° first-order derivative filters More complicated transforms can also be incorporated in the proposed framework in a straightforward manner (possibly at the expense of computational complexity).

Fig. 1

Contourlets in (a) horizontal and (b) vertical directions.

Using these transforms, the image model can be constructed to exploit the sparsity in the transform coefficients. A commonly used sparse image model is [20] where denotes the l-pseudonorm, and α are the weighting coefficients. It is known from the compressive sensing and sparse representation literature [15, 16, 20] that using 0 < p ≤ 1 enforces sparsity on Df, with smaller p values increasing the sparsity effect. Hence, the prior in Eq. (12) enforces sparsity in the transform coefficients Df, which in turn leads to smoothness in the image estimate.

The disadvantage of using the prior in Eq. (12) is that it is nonconvex, and using this penalty creates a high number of local minima when estimating f. Instead, we use separate Gaussian priors on each transform coefficient or in a more compact form as where A are diagonal matrices with α, i = 1,... N in the diagonal. Compared to Eq. (12), where a single parameter is assigned to all coefficients of kth filter output, separate parameters are used for each coefficient. It can be shown that Eq. (14) is equivalent to Eq. (12) in the limit p → 0 [21], hence it highly enforces sparsity. The model Eq. (14) has the advantage of being convex (as opposed to Eq. (12)), and therefore optimization over Eq. (14) is much easier and more robust compared to l minimization.

The parameters α have a special important role in Eq. (13): they represent the local spatial activity at each location, and hence they are a measure of spatial variation in the corresponding filters direction. It is clear that the model Eq. (14) requires a large number of parameters, whose manual selection is not practical. We can, however, estimate them simultaneously with the complex image. For their estimation, we employ an additional level of model and assign uniform priors Notice that this modeling assigns equal probability to all possible values of α, hence no prior knowledge is assumed on its value.

It should be emphasized that this modeling based on sparsity principles does not necessarily cause the estimates to have very sparse coefficients. Real images are generally only approximately sparse, i.e., a few transform coefficients have large values whereas most coefficients are very small, but not necessarily exactly zero. This behavior is generally referred to as compressible [22]. Enforcing sparsity to an extreme extent can therefore suppress subtle image features, which may be important. Our modeling in Eq. (14), on the other hand, allows for the small-valued transform coefficients through Gaussian distributions while enforcing the general compressible structure of the images.

4.2. Noise model

The signal-independent noise is modeled via a independent Gaussian noise model on the observed field as with σ2 the noise variance. An additional level of modeling (as in Eq. (15)) can be incorporated to estimate this parameter as well. However, SLIM provides images with very high SNRs (on the order of 1000 or more), and therefore σ2 is generally very small and can be estimated experimentally from an uniform area of the observed image. In addition, the Gaussian noise assumption becomes an accurate description of the noise in SLIM due to the high SNR (as in fluorescence microscopy [10]).

4.3. Algorithm

Using the modeling described in the previous sections, we formulate the problem of estimating the unknown complex image f and the parameters α using the maximum a posteriori (MAP) estimates, that is, This problem is convex in f and α, but not jointly, and therefore we resort to an iterative scheme to estimate the unknowns in an alternating fashion. The optimal estimate of the complex image can be found by taking the derivative of Eq. (18) and setting it equal to zero, which results in

The parameters α can be estimated in a similar way by equating the corresponding derivatives to zero, which results in where ε is a small number (e.g., 10−6) used to avoid numerical instability. It follows from Eq. (20) that the parameters α are functions of the kth filter response at pixel i, and therefore a spatially-adaptive estimation is employed for f in Eq. (19) through their joint estimation. Notice also that matrices A are spatial-adaptivity matrices controlling the smoothness applied at each location; when the filter responses at a pixel are very small, the algorithm assumes that the pixel has low spatial variation, and applies a large amount of smoothness at that point. On the other hand, if the filter responses are high, the pixel is likely to be close to an edge and the smoothness amount is lowered to preserve the image structure.

In summary, the proposed method dSLIM consists of alternating estimations of the complex field f using Eq. (19), and the spatial adaptivity matrices A using Eq. (20). The block diagram of a single dSLIM iteration is shown in Fig. 2. The filters D consist of the derivative operators in Eq. (10) and Eq. (11), and the directional contourlets shown in Fig. 1. The estimate of the image Eq. (19) can be computed very efficiently using the conjugate gradient (CG) method. The matrices H and D do not have to be explicitly constructed during CG iterations; all operations in Eq. (19) can be performed via convolutions in the spatial domain or multiplications in the Fourier domain. Empirically, we found that the algorithm converges rapidly; a few iterations (up to 5–10) is generally enough to provide high-quality results. Hence, the proposed method can be applied to large images very efficiently.

Fig. 2

Block diagram of one dSLIM iteration.

The proposed algorithm contains only one free parameter (the noise variance σ2) that needs to be set by the user. In our experiments, we empirically estimated its value by taking a rectangular region of the image with uniform values and computing the variance in this region. As mentioned above, this estimate is known to be reliable in images with high SNR [10]. Since this free parameter corresponds to a physical quantity, its estimation is relatively easy and does not require extensive image-dependent tuning.

Finally, it should be emphasized that dSLIM does not alter the quantitative imaging property of SLIM, which is one its main advantages. As the deconvolution is applied to the complex field exp [iΦ̃(x,y)] rather than to the measured intensities I(x,y,ϕ0) in Eq. (3), the quantitative phase information is preserved. In contrast, traditional deconvolution methods [10] applied directly to the intensity images can not preserve the quantitative information.

5. Experiments

In this section, we illustrate the application of dSLIM to complex field images obtained by SLIM and quantitatively demonstrate the resolution increase. All SLIM images were acquired using a white-light source (mean wavelength λ = 530 nm); the field of view is 75μm × 100μm with the CCD resolution of 1040 × 1388. In all reported experiments, the specimen is relatively thin such that the whole image is in focus, and the degradation in the image is only due to a planar PSF. The PSF, depicted in Fig. 3(a), is obtained experimentally by imaging a sub-resolution 200nm microbead treated as a point-source. Due to the high SNR provided by SLIM, this PSF closely matches the actual optical transfer function of the imager.

Fig. 3

(a) Experimental PSF, (b) result of dSLIM, (c) normalized horizontal cross-sections, and (d) estimated FWHMs during the iterative procedure.

In all images, the noise level is estimated within the range 10−7–10−6 (for a maximum signal value of 1), which is used as the value of the parameter σ2. The NAs of the objective and condenser are NAo = 0.75 and NAc = 0.55, respectively. The experimentally measured full-width-at-half-maximum (FWHM) of the PSF is 540 nm, which is comparable with the expected Rayleigh limit, calculated as .

We first investigate the resolution increase obtained by dSLIM by applying it to the experimental PSF. The experimental PSF is shown in Fig. 3(a). Treating this as the original image, we apply dSLIM and obtain the result shown in Fig. 3(b). The FWHM of the original PSF is approximately 540nm, whereas after dSLIM, the FWHM is reduced to approximately 238nm, corresponding to a 2.3 times increase in resolution. The horizontal cross-section of the images are shown in Fig. 3(c), where the reduction in FWHM is clearly visible. Notice that this resolution is significantly below the diffraction limit. The estimated FWHM during the iterative procedure is shown in Fig. 3(d), which shows a significant reduction in the first iterations and convergence within 10 iterations.

To further examine the resolution increase provided by dSLIM, we next apply it to SLIM images of multiple 200nm microbeads. Figure 4(a) shows an image of two beads approximately 550 nm apart, which are barely resolved in the original SLIM image. The image after applying dSLIM is shown in Fig. 4(b), where the microbeads are clearly separated while their distance is accurately preserved (Fig. 4(c)).

Fig. 4

Images of two microbeads (a) SLIM and (b) dSLIM. The cross-sectional profiles (with normalized maximum phase values) are shown in (c).

Next, we demonstrate dSLIM images of biological specimen. A SLIM phase image of a hippocampal neuron is shown in Fig. 5. The SLIM images are shown on the left column, while dSLIM images are shown on the right column. It is clear that dSLIM effectively removes the blur, while deconvolution artifacts and noise are successfully suppressed and object boundaries faithfully preserved. dSLIM recovers the details of the fine structure of the specimen which are hard to observe in the original image. The high quality of the dSLIM result can be better observed from the detailed parts of the image shown in Fig. 5(c) in comparison with the SLIM image.

Fig. 5

(a) SLIM image of a hippocampal neuron, (b) image provided by dSLIM. Detailed images of the central parts are shown in (c) and (d).

Another example is shown in Fig. 6 ( Media 1, Media 2), which is a snapshot of a dynamic SLIM image sequence of a live hippocampal neuron culture. Interference-microscopy is extremely useful in monitoring dynamic cellular processes over time, as it does not require invasive contrast enhancement techniques (such as fluorescence tagging). SLIM is a very attractive modality for this application due to its very high spatial and temporal resolution (several frames/second). It is clear from Fig. 6(c) (left) ( Media 2) that the observation noise level is very low, but the image exhibits a certain level of blur due to the diffraction-limited PSF. The result of dSLIM is shown in Fig. 6(c) (right) ( Media 2), which shows a clear resolution improvement over the original image. The increased spatial resolution also positively affects the examination of dynamic neuron processes. Due to more accurate estimates of size and locations of the particles, the dynamic changes and hence the biological behavior can be better observed (see Media 2 for a visualization).

Fig. 6

SLIM dynamic imaging of live hippocampal neuron in primary cell culture. (a) SLIM image ( Media 1), (b) dSLIM image ( Media 1), (c) Detailed areas of regions A, B, C from (a) and (b) ( Media 2). The SLIM image regions are shown on the left, while the dSLIM image regions are on the right.

To further confirm the increase in resolution in real images, we examine microparticles in the SLIM image shown in Fig. 6. Figures 7(a) and 7(b) show detailed images of a single particle from the SLIM and dSLIM images (marked as region D in Fig. 6(a)). The cross-sections of the images are shown in Fig. 7(c). The vertical diameter of the particle is measured as approximately 1.5 μm in the SLIM image, whereas it is measured as 0.63 μm in the dSLIM image. The reduction in the length of the particle is approximately 2.3, which is in agreement with the result of the PSF deconvolution experiment (Fig. 3).

Fig. 7

(a) SLIM image of a single particle from the region D in Fig. 6, (b) dSLIM image, and (c) normalized cross-sections from the images.

Our final experiment shows two neuronal processes (putative axons) which were not resolved in the original SLIM image (region marked as E in Fig. 6(a)). The detailed area is shown in Fig. 8(a), and the dSLIM result is shown in Fig. 8(b). dSLIM reveals two objects located approximately 430 nm apart. This can also be observed from the normalized cross-sections passing through the maximum phase values, shown in Figs. 8(c). The evolution of this area over time is shown in Fig. 9 ( Media 3). It can be observed both from the original images and the cross-sections that the objects are just resolved in some time frames, but unresolved in others. dSLIM successfully separates the objects through the whole dynamic sequence ( Media 3).

Fig. 8

(a) Two very closely located particles from the region E in Fig. 6 not resolved in the SLIM image, (b) dSLIM image, and (c) the normalized cross-sections through the maximum phase values.

Fig. 9

Dynamic evolution of the area shown in Fig. 8 ( Media 3). Top row: SLIM images, middle row: dSLIM images, bottom row: Normalized cross-sections of the images (through the segment shown in the top-left image) at each time point.

6. Conclusion

In this paper, we presented a novel deconvolution method, dSLIM, for complex image fields acquired by interference microscopy. Our formulation is based on three key observations: First, the image formation can be treated as a linear process in the complex fields, such that the degradation of the microscopy can be modeled by a PSF acting on the complex images. Second, due to the high SNR provided by the SLIM, the experimentally obtained PSF of the imager closely matches the actual PSF. Finally, the phase images of biological specimen can be very accurately modeled using sparsity principles. We combined these properties to develop a very effective deconvolution procedure that significantly improves the final resolution, allowing imaging very fine structures and motions in live cells below the diffraction limit. Due to the high spatial and temporal resolution, this approach can be utilized to acquire new information for studying live cells.