Light field microscopy (LFM) is a scanning-free 3-dimensional (3D) microscopic imaging method based on the theory of light field. This technique allows sub-second (~10 Hz) large volumetric imaging ([~0.1 to 1 mm]3) with ~1 μm spatial resolution in the condition of weak scattering and semi-transparence, which has never been achieved by other methods. Just as in traditional light field rendering, there are two steps for LFM imaging: light field capture and processing. In most setups, a microlens array is used to capture the light field. As for processing, it can be based on two kinds of representations of light propagation: the ray optics picture and the wave optics picture. The Stanford University Computer Graphics Laboratory published their first prototype LFM in 2006 and has been working on the cutting edge since then.
Light field generation.
A light field is a collection of all the rays flowing through some free space, where each ray can be parameterized with four variables. In many cases, two 2D coordinates–denoted as formula_1 &amp; formula_2–on two parallel planes with which the rays intersect are applied for parameterization. Accordingly, the intensity of the 4D light field can be described as a scalar function: formula_3, where formula_4 is the distance between two planes.
LFM can be built upon the traditional setup of a wide-field fluorescence microscope and a standard CCD camera or sCMOS. A light field is generated by placing a microlens array at the intermediate image plane of the objective (or the rear focal plane of an optional relay lens) and is further captured by placing the camera sensor at the rear focal plane of the microlenses. As a result, the coordinates of the microlenses formula_1 conjugate with those on the object plane (if additional relay lenses are added, then on the front focal plane of the objective) formula_6; the coordinates of the pixels behind each microlens formula_2 conjugate with those on the objective plane formula_8. For uniformity and convenience, we shall call the plane formula_6 the original focus plane in this article. Correspondingly, formula_4 is the focal length of the microlenses (i.e., the distance between microlens array plane and the sensor plane).
In addition, the apertures and the focal-length of each lens and the dimensions of the sensor and microlens array should all be properly chosen to ensure that there is neither overlap nor empty areas between adjacent subimages behind the corresponding microlenses.
Realization from the ray optics picture.
This section mainly introduces the work of Levoy "et al"., 2006.
Perspective views from varied angles.
Owing to the conjugated relationships as mentioned above, any certain pixel formula_11 behind a certain microlens formula_12 corresponds to the ray passing through the point formula_13 towards the direction formula_14. Therefore, by extracting the pixel formula_11 from all subimages and stitching them together, a perspective view from the certain angle is obtained: formula_16. In this scenario, spatial resolution is determined by the number of microlenses; angular resolution is determined by the number of pixels behind each microlens.
Tomographic views based on synthetic refocusing.
Step 1: Digital refocusing.
Synthetic focusing uses the captured light field to compute the photograph focusing at any arbitrary section. By simply summing all the pixels in each subimage behind the microlens (equivalent to collecting all radiation coming from different angles that falls on the same position), the image is focused exactly on the plane that conjugates with the microlens array plane:
formula_17,
where formula_18 is the angle between the ray and the normal of the sensor plane, and formula_19 if the origin of the coordinate system of each subimage is located on the principal optic axis of the corresponding microlens. Now, a new function can defined to absorb the effective projection factor formula_20 into the light field intensity formula_21 and obtain the actual radiance collection of each pixel: formula_22.
In order to focus on some other plane besides the front focal plane of the objective, say, the plane whose conjugated plane is formula_23 away from the sensor plane, the conjugated plane can be moved from formula_4 to formula_25 and reparameterize its light field back to the original one at formula_4:
formula_27.
Thereby, the refocused photograph can be computed with the following formula:
formula_28.
Consequently, a focal stack is generated to recapitulate the instant 3D imaging of the object space. Furthermore, tilted or even curved focal planes are also synthetically possible. In addition, any reconstructed 2D image focused at an arbitrary depth corresponds to a 2D slice of a 4D light field in the Fourier domain, where the algorithm complexity can be reduced from formula_29 to formula_30.
Step 2: Point spread function measurement.
Due to diffraction and defocus, however, the focal stack formula_31 differs from the actual intensity distribution of voxels formula_32, which is really desired. Instead, formula_31 is a convolution of formula_32 and a point spread function (PSF):
 formula_35
Thus, the 3D shape of the PSF has to be measured in order to subtract its effect and to obtain voxels' net intensity. This measurement can be easily done by placing a fluorescent bead at the center of the original focus plane and recording its light field, based on which the PSF's 3D shape is ascertained by synthetically focusing on varied depth. Given that the PSF is acquired with the same LFM setup and digital refocusing procedure as the focal stack, this measurement correctly reflects the angular range of rays captured by the objective (including any falloff in intensity); therefore, this synthetic PSF is actually free of noise and aberrations. The shape of the PSF can be considered identical everywhere within our desired field of view (FOV); hence, multiple measurements can be avoided.
Step 3: 3D deconvolution.
In the Fourier domain, the actual intensity of voxels has a very simple relation with the focal stack and the PSF:
formula_36,
where formula_37 is the operator of the Fourier transform. However, it may not be possible to directly solve the equation above, given the fact that the aperture is of limited size, resulting in the PSF being bandlimited (i.e., its Fourier transform has zeros). Instead, an iterative algorithm called "constrained iterative deconvolution" in the spatial domain is much more practical here: 
This idea is based on constrained gradient descent: the estimation of formula_32 is improved iteratively by calculating the difference between the actual focal stack formula_31 and the estimated focal stack formula_42 and correcting formula_32 with the current difference (formula_32 is constrained to be non-negative).
Fourier Slice Photography.
The formula of formula_45 can be rewritten by adopting the concept of the Fourier Projection-Slice Theorem. Because the photography operator formula_46 can be viewed as a shear followed by projection, the result should be proportional to a dilated 2D slice of the 4D Fourier transform of a light field. Precisely, a refocused image can be generated from the 4D Fourier spectrum of a light field by extracting an 2D slice, applying an inverse 2D transform, and scaling. Before the proof, we first introduce some operators:
By these definitions, we can rewrite formula_53.
According to the generalized Fourier-slice theorem, we have
formula_54,
and hence the photography operator has the form
formula_55.
According to the formula, we know a photograph is the inverse 2D Fourier transform of a dilated 2D slice in the 4D Fourier transform of the light field.
Discrete Fourier Slice Photography.
If all we have available are samples of the light field, instead of use Fourier slice theorem for continuous signal mentioned above, we adopt discrete Fourier slice theorem, which is a generalization of the discrete Radon transform, to compute refocused image.
Assume that a lightfield formula_56 is periodic with periods formula_57 and is defined on the hypercube formula_58. Also, assume there are formula_59 known samples of the light field formula_60, where formula_61 and formula_62, respectively. Then, we can define formula_63 using trigonometric interpolation with these sample points:
formula_64,
where
formula_65.
Note that the constant factors are dropped for simplicity.
To compute its refocused photograph, we replace infinite integral in the formula of formula_66 with summation whose bounds are formula_67 and formula_68. That is,
formula_69.
Then, by discrete Fourier slice theorem indicates, we can represent the photograph using Fourier slice:
formula_70
Realization from the wave optics picture.
Although ray-optics based plenoptic camera has demonstrated favorable performance in the macroscopic world, diffraction places a limit on the LFM reconstruction when staying with ray-optics parlance. Hence, it may be much more convenient to switch to wave optics. (This section mainly introduce the work of Broxton "et al"., 2013.)
Discretization of the space.
The interested FOV is segmented into formula_71 voxels, each with a label formula_72. Thus, the whole FOV can be discretely represented with a vector formula_73 with a dimension of formula_74. Similarly, a formula_75 vector formula_76 represents the sensor plane, where each element formula_77 denotes one sensor pixel. Under the condition of incoherent propagation among different voxels, the light field transmission from the object space to the sensor can be linearly linked by a formula_78 measurement matrix, in which the information of PSF is incorporated:
 formula_79
In the ray-optics scenario, a focal stack is generated via synthetically focusing of rays, and then deconvolution with a synthesized PSF is applied to diminish the blurring caused by the wave nature of light. In the wave optics picture, on the other hand, the measurement matrix formula_80–describing light field transmission–is directly calculated based on propagation of waves. Unlike transitional optical microscopes whose PSF shape is invariant (e.g., Airy Pattern) with respect to position of the emitter, an emitter in each voxel generates a unique pattern on the sensor of a LFM. In other words, each column in formula_80 is distinct. In the following sections, the calculation of the whole measurement matrix would be discussed in detail.
Optical impulse response.
The optical impulse response formula_82 is the intensity of an electric field at a 2D position formula_83 on the sensor plane when an isotropic point source of unit amplitude is placed at some 3D position formula_84 in the FOV. There are three steps along the electric-field propagation: traveling from a point source to the native image plane (i.e., the microlens array plane), passing through the microlens array, and propagating onto the sensor plane.
Step 1: Propagation cross an objective.
For an objective with a circular aperture, the wavefront at the native image plane formula_85 initiated from an emitter at formula_86 can be computed using the scalar Debye theory:
formula_87,
where formula_88 is the focal length of the objective; formula_89 is its magnification. formula_90 is the wavelength. formula_91 is the half-angle of the numerical aperture (formula_92 is the index of refraction of the sample). formula_93 is the apodization function of the microscope (formula_94 for Abbe-sine corrected objectives). formula_95 is the zeroth order Bessel function of the first kind. formula_96 and formula_97 are the normalized radial and axial optical coordinates, respectively:
formula_98
formula_99,
where formula_100 is the wave number.
Step 2: Focusing through the microlens array.
Each microlens can be regarded as a phase mask:
formula_101,
where formula_102 is the focal length of microlenses and formula_103 is the vector pointing from the center of the microlens to a point formula_104 on the microlens. It is worth noticing that formula_105 is non-zero only when formula_104 is located at the effective transmission area of a microlens.
Thereby, the transmission function of the overall microlens array can be represented as formula_107 convoluted with a 2D comb function:
formula_108,
where formula_109 is the pitch (say, the dimension) of microlenses.
Step 3: Near-field propagation to the sensor.
The propagation of wave front with distance formula_4 from the native image plane to the sensor plane can be computed with a Fresnel diffraction integral:
formula_111,
where formula_112 is the wave front immediately passing the native imaging plane.
Therefore, the whole optical impulse response can be expressed in terms of a convolution:
formula_113.
Computing the measurement matrix.
Having acquired the optical impulse response, any element formula_114 in the measurement matrix formula_115 can be calculated as:
formula_116,
where formula_117 is the area for pixel formula_118 and formula_119 is the volume for voxel formula_120. The weight filter formula_121 is added to match the fact that a PSF contributes more at the center of a voxel than at the edges. The linear superposition integral is based on the assumption that fluorophores in each infinitesimal volume formula_122 experience an incoherent, stochastic emission process, considering their rapid, random fluctuations.
Solving the inverse problem.
The noisy nature of the measurements.
Again, due to the limited bandwidth, the photon shot noise, and the huge matrix dimension, it is impossible to directly solve the inverse problem as: formula_123. Instead, a stochastic relation between a discrete light field and FOV more resembles:
formula_124,
where formula_125 is the background fluorescence measured prior to imaging; formula_126 is the Poisson noise. Therefore, formula_127 now becomes a random vector with Possion-distributed values in units of photoelectrons e−.
Maximum likelihood estimation.
Based on