1. Introduction

Conventional cameras or microscopes capture two-dimensional (2D) image of the scene, whereas the depth (z) information is missing as the acquisition essentially integrates the light field along the z direction. However, the depth information is very important in many applications such as entertainment, optical inspection, biological imaging, and so on. In the last decade, a rapid, scanless three-dimensional (3D) imaging approach, called light field imaging, has received more and more attention because of its high compatibility with various applications ranging from photographic camera [1] to microscope [2, 3]. From the geometric optics point of view, light field imaging can simultaneously record the 2D spatial and 2D angular information of each light ray in the bundle of rays observed by the digital image sensors, thus allowing the generation of perspective views and focal stacks of the scene by computationally tracing the rays back to certain planes. Technically, the light field can be captured by using either a microlens array with a standard camera [1, 2], a camera array [4, 5], or an absorbing mask [6]. Of all these, microlens-array-based light field instruments receive particular interest because of its commercial success and potential application in light field microscopy [2, 3, 7]. However, there is an intrinsic trade-off between the spatial and the angular resolution due to the microlens array: one must sacrifices spatial resolution in order to obtain angular resolution, which determines the number of perspective views, number of focal stacks, and the accuracy of the estimation of depth map.

Alternatively, one can manage to collect the angular moment of the light rays using the light field moment imaging (LMI) technique recently proposed by Orth and Crozier [8]. With the LMI, one does not need to sacrifices the spatial resolution to achieve the perspective views since the raw images are taken by standard cameras or light microscopes. Recently, the LMI has also been used in Fourier hologram synthesis [9] and coherent diffractive imaging [10]. The principle of the LMI is to extract the first-order angular moment of the light field by solving the Poisson equation, which connects the Laplacian of the scalar potential of the light field with the intensity derivative, in a way similar to the transport of intensity equation [11]. We note that the Poisson solver severely affected by how accurately the intensity derivative can be estimated. To obtain a good estimation, one can think of taking the two intensity images at two lateral plane as close as possible from the pure mathematical perspective. However, as the two planes are getting close to each other, noise of the sensors becomes increasingly significant, and eventually submerges the intensity derivative, leading to a fault reconstruction.

In this paper, we propose a light field moment microscopy (LFMM) inspired by the work of Orth and Crozier [8] for the inspection of samples in micro-scale. Significant difference between our LFMM and the standard LMI is that we use a better method to estimate the intensity derivative in the Poisson solver. Both numerical and optical experiments are carried out to demonstrate our method. The paper is organized as follows: in Sec. 2, we review the basic principle of LMI. Then, we discuss the noise problem in LMI and describe our method of noise reduction LMI in Sec. 3. The simulation and experiment results of the proposed approach is described in Sec. 4 and Sec. 5, respectively. Finally, the conclusion is drew in Sec 6.

2. Light field moment imaging

According to the plenoptic function [12], the light field can be parameterized as a five-dimensional function L(x, y, z, tanθ_X, tanθ_Y), where (x, y, z) is the spatial coordinates and (tanθ_X, tanθ_Y) is the angular coordinates. The function L describes the position and propagation direction of each ray in the light field. The lateral intensity I(x, y) observed at a plane z can be expressed as the integral of L with respect to the angular dimensions [13] and over the distance z. The LMI connects the derivative of I(x, y; z) with respect to z to a scalar potential U(x, y; z) as [8]:

It is clear that in the coherent limit, Eq. (1) reduces to the transport of intensity equation, and the scalar potential U becomes the phase of the field. Numerically, Eq. (1) can be solved using a Fourier-tranform-based algorithm

Once the scalar potential U(x, y; z) is computed [Eq. (3)], we can reconstruct the angular moment vector field M via the definition of U (x, y; z)

Orth and Crozier have shown that σ ² = NA ² gives a good perspective shifting effect [8]. By extracting the 2D slice of the light field L at different (tanθ_X, tanθ_Y) coordinate, we have different perspective view of the scene, similar to the perspective-shifting effect observed in the light-field photography [1] and light-field microscopy [2].

3. Noise reduction LFMM

It is seen from Eq. (5) in Section 2 that the intensity derivative ∂ I/∂ z is estimated from the intensity measured at the two neighoring planes. Ideally, the estimation accuracy is determined by the distance Δz between the two planes according to the finite difference theory [14]. Thus, Δz should be as small as possible in order to obtain a good estimation. Otherwise, the influence of the higher-order terms cannot be ignored, and the linearity assumption to finite difference approximation is invalid, leading to a significant noise in the reconstructed light field. Therefore, the most direct way to increase the estimation accuracy of ∂ I/∂ z is to capture the two images at two planes as close as possible. However, it becomes difficult to distinguish the difference between I ₁ (x, y; z ₁) and I ₂ (x, y; z ₂) when the two planes are too close owing to the dynamic range and noise characteristics of the camera in actual experiments, in particular when the object is weakly diffractive. As a result, significant noise appears in the difference I ₁ (x, y; z ₁) − I ₂ (x, y; z ₂). Furthermore, the Poisson solver described in Eq. (3) requires the construction of the function H(f_x, f_y) defined in Eq. (4). This filter function has an extremely large value as $f_x^2+f_y^2\to 0$, and is therefore very sensitive to the noise at lower spatial frequency region. This means that the reconstructed scalar potential U (x, y; z) is vulnerable to noise at this area as well, limiting the effect of the perspective view in LMI.

It becomes clear now that the approach in [8] Orth and Crozier proposed must have a trade-off between the estimation accuracy of ∂ I/∂ z and noise. In this work, we propose to fit the intensity derivative using multi-images to address this problem. In this algorithm, the intensity derivative ∂ I/∂ z is calculated in a pixelwise manner. At each pixel (x_i,y_j), where i, j = 1,…N, we first fit the intensity evolution to an m ^th order polynomial. Then the intensity derivative at this pixel, ∂ I(x_i, y_j)/∂ z, is calculated by directly taking the first order coefficient of the polynomial. One can see that all the images contribute to the estimation of intensity derivative, which will help improve the estimation accuracy. The similar method has been proposed in the field of phase retrieval for improving the accuracy of phase reconstruction [15].

Now let us provide a more theoretical analysis to the problem. To begin, we take the Taylor series expansion of the intensity captured at a defocus distance Δz, and have [14]

Apparently, we need to estimate the higher-order terms in order to improve the estimation accuracy of ∂ I/∂ z. It has been shown that the n^th order derivative of a function can be computed numerically from n + 1 or more equally spaced measurements using polynomial fitting. Then we can use the more accurate estimation of ∂ I/∂ z in the Poisson solver, and Eq. (3) now can be written as

It is expected to obtain a better solution by solving this higher order Poisson equation. In the following sections, we provide both numerical simulation and experimental results to demonstrate this method.

4. Numerical simulation

In numerical simulation, we used a pure phase object to verify the usefulness of the higher order terms as this will provides a more quantitative understanding of the method. To put it simple, we use a coherent imaging system in our simulation. As discussed above, in the coherent limit, the light field moment reduces to the spatial frequency. In simulations, we first numerically construct a 4f imaging system based on wave optics, and propagate the object wavefront to a number of axial planes at both sides of the focal plane by using the Fresnel transform [16, 17], and store the intensity obtained at each of these planes. The distance between every two neighboring planes is set to the same value. Using this intensity data set, we can examine the evolution of the intensity at each pixel of these planes, and fit it to a curve specified by a high-order polynomial. The reduced noise ∂ I/∂ z can be calculated by directly taking the first order coefficient of the polynomial.

In simulation, we obtained intensity images with the axial spacing of Δz = 1 μm by using our numerical coherent imaging system. The simulation results are shown in Fig. 1 and Fig. 2. Figure 1(a) shows the axial evolution of the intensity of a certain pixel, fitted to the 1st, the 7th and the 13th orders polynomial, respectively. This clearly shows that the estimation of the intensity derivative using the first-order difference [I(Δz) − I(0)]/Δz only gives a very rough estimation, with significant deviation from the true value (shown by the blue markers). Instead, high-order polynomial fitting gives a far better estimation. To quantify the accuracy, we plot in Fig. 1(b) the fit error, the value of which corresponding to the 1st, 7th and 13th order polynomial fitting is 0.0249, 0.0065 and 0.0008 RMS, respectively. It is therefore expected to recover the moment M and the light field L using the better estimation of the intensity derivative, i.e., solving Eq. (10) instead of Eq. (3).

 

Fig. 1 Simulation results: (a) Change of the intensity of the field at a certain pixel due to the propagation along the z direction (blue markers), and the fitting polynomial to the 1st, the 7th and the 13th order, and (b) the corresponding fit errors.

Download Full Size | PDF

 

Fig. 2 The s (first row) and t (second row) component of the moment vector field M. The two images in the first column show the ground truth data for comparison, and the images in the other three columns are recovered using the intensity derivative estimated by fitting intensity variation along z to the 1st, the 7th and the 13th orders polynomials, respectively.

Download Full Size | PDF

Poisson noise was added to the simulated intensity patterns in order to test the noise reduction performance of our LFMM technique. The recovered moment vectors M using different orders of polynomial fitting are plotted in Fig. 2. The first column in Fig. 2 shows the s and t components of the ideal moment M. These two images serve as the ground-truth data to test the proposed method. The images in the second column were computed by solving Eq. (3). And the images in the third and the fourth columns were computed by solving the high-order LMI equation [Eq. (10)] with the intensity derivative fitted to the 7th and the 13th order, respectively.

It is clearly seen from Fig. 2 that the estimation of the intensity derivation using the finite difference gives the worst result due to the noise sensitivity of the LMI solver [8]. The result is significantly improved when using the higher-order estimation of the intensity derivative. However, we observed in the third and the fourth columns that the angular moments recovered using the 13th order polynomial fitting is worse than the one recovered using the 7th order fitting. The reason for this is that the higher-order polynomial fitting may fit the noise as the signal, yielding smaller RMS errors. This is the overfitting issue [18] encountered in curve fitting problems.

5. Experimental results

Now we set up an experimental system as shown in Fig. 3(a) to demonstrate the proposed LFMM technique. The system is just a standard microscope (Nikon Ni) mounted with a 10 × 0.3 NA microscope objective (Nikon) and an sCMOS camera (PCO Edge 4.2 M). The sample we used was a mosquitoes’ mouth-parts, which was a common biological sample slice one can buy on the market. We captured a stack of intensity images of the sample with an axial spacing of Δz = 1 μm by tuning the focus knobs. Three of intensity images are shown in Fig. 3(b).

 

Fig. 3 (a) Experimental setup. The intensity focal stack images of the mosquitoes’ mouth-parts on a microscopic slice are captured by adjusting the focus knob. (b) Three of the sample images that we captured.

Download Full Size | PDF

It is clearly seen from Fig 3(b) that the intensity images captured at adjacent planes are very close to each other. This means that the corresponding intensity difference is easily corrupted by noise, as evidenced in Fig. 4(a). Instead, the noise can be significantly suppressed by using the proposed method, as shown in Fig. 4(b)–Fig. 4(d), in which the intensity derivative was calculated by fitting the polynomial to the 5th, the 7th and the 13th order using the experimental data set.

 

Fig. 4 Intensity derivative ∂ I/∂ z calculated by using (a) the finite difference [Eq. (5)], and (b) the 5th, (c) the 7th, and (d) the 13th order fitting polynomials [Eq. (9)].

Download Full Size | PDF

With these data at hand, we can recover the scalar potential U and the moment M by solving Eq. (10) and Eq. (6). The experimental results are plotted in Fig. 5. We observed that the recovered potentials and moments contain more detailed and finer structures by using our method (the second and third column) in comparison to the results recovered using the finite difference (the first column). This suggests that the proposed LFMM can significantly reduce the noise. In our experiments, we also observed the effect of overfitting in the recovered U and M. This is shown in the fourth column in Fig. 5 as well as in Fig. 4(d). One can see that noise is present in these images because fitting the data to a polynomial of too high order will take the noise as signal.

 

Fig. 5 Experimental results of the recovered scalar potential U (first row) and the s (second row) and t (third row) component of the moment vector field M using the intensity derivative estimated by fitting the intensity variation along z to different order of polynomial.

Download Full Size | PDF

We can then construct the light field L of the scene according to Eq. (7) after recovering the moment vector M. Different perspective views of the sample can then be achieved by extracting the corresponding 2D intensity slice out of the constructed light field L. Examples of the extracted perspective views of the mosquitoes’ mouth-parts at angles coordinates (0°, 5.8°) are shown in Fig. 6. One can find in Visualization 1, Visualization 2, Visualization 3, and Visualization 4 for the complete perspective view of the light field. The constructed light-fields shown in Fig. 6 and the Visualizations suggest that the proposed LFMM technique has a much better overall noise performance in comparison to the standard LMI. The noise performance of the LFMM is dependent on the estimation accuracy of the intensity derivative. In our experiments, we observed that the intensity derivative estimated using the 7th-order fitting offers the best perspective-shifting impression (see Visualization 3). The other two are either suffered from noise (see Visualization 2) due to underfitting, or suffered from degeneration, especially on the marginal perspective view (see Visualization 4), by virtue of overfitting.

 

Fig. 6 2D intensity slice extracted from the constructed 4D light field with the angular coordinates (θ_X, θ_Y) equal to (0°, 5.8°). It is clearly seen that the 7th-order fitting offers the best performance. See Visualization 1, Visualization 2, Visualization 3, and Visualization 4 for the view of the perspective shifting effect.

Download Full Size | PDF

6. Conclusion

To summarize, we have described a noise reduction LFMM technique. The proposed technique makes a better estimation of the intensity derivative in solving the LMI equation by using polynomial fitting instead of finite difference, resulting in the significant reduction of noise. We have performed both simulation and experiments to demonstrate the proposed technique. From the experimental results, we can clearly see the light-field-like perspective view of the biological sample by using a standard microscope, without the sacrifice of the spatial resolution as in standard light field microscopy [2,3,7]. We believe that the proposed technique will open up new possibility for standard microscopes in high-resolution light-field observations.

It is also necessary to discuss the computational complexity of the noise reduction algorithm. Note that the computational complexity of polynomial fitting using the least square method is O(M ² P), where M is the number of points, and P is the order of polynomial to be fit. So, for the images of N × N pixels as in our case, the total computation cost is O(M ² N ² P). It is worthy of pointing out that the pixelwise treatment of the algorithm can be implemented using parallel computing. In that case, the computation can be significantly sped up. In a follow-up study, we will try to find a better approach to reduce the noise with fewer images and lighter computation load.