Artifact-free deconvolution in light field microscopy

: The sampling patterns of the light ﬁeld microscope (LFM) are highly depth-dependent, which implies non-uniform recoverable lateral resolution across depth. Moreover, reconstructions using state-of-the-art approaches suﬀer from strong artifacts at axial ranges, where the LFM samples the light ﬁeld at a coarse rate. In this work, we analyze the sampling patterns of the LFM, and introduce a ﬂexible light ﬁeld point spread function model (LFPSF) to cope with arbitrary LFM designs. We then propose a novel aliasing-aware deconvolution scheme to address the sampling artifacts. We demonstrate the high potential of the proposed method on real experimental data.


Introduction
Since it was first proposed by Levoy et al. in 2006 [1], light field microscopy has proven very useful in biological applications involving fast dynamics, due to its high speed 3D imaging capability. The potential of the modality has been show-cased in various applications, including live cell imaging [2], fish larvae neuro-dynamics recording [3,4], heart imaging and blood flow monitoring [5].
The light field microscope (LFM) enables scan-less 3D imaging of fluorescent specimens by incorporating an array of micro-lenses into the optical path of a conventional wide-field microscope. Thus, both spatial and directional light field information is captured in a single shot, allowing for subsequent volumetric reconstruction of the imaged sample.
Based on the integral imaging principles [6], the lenslet-based (plenoptic) devices [7] have established plenoptic imaging as an active topic in computational imaging, enabling postacquisition refocusing [8][9][10] or 3D imaging of transparent microscopic samples [1,11]. In the early stages, the methods for rendering extended depth of field images from plenoptic devices were limited to lenslet resolution [8], which is the number of available micro-lenses. In recent years, various approaches to address this limitation were proposed. Examples include variations of the hardware to allow for attenuation masks [12] or wavefront coding techniques [13], which demonstrated higher recoverable resolutions and extended depth of field. On the other hand, inspired by the large amount of work on computational super-resolution in the computer vision field [14,15], algorithms for super-resolving the light field were developed, involving multi-view reconstruction [16,17], or explicit image formation models for plenoptic devices employing either ray-based [18][19][20] or wave-based optics [13,21,22].
In [22], Broxton et al. introduced a wave-based model to describe the propagation of light through an original plenoptic LFM setup [1], together with a 3D deconvolution method. They demonstrate superior reconstructions in terms of lateral resolution (compared to lenslet resolution) for most of the axial range. The improvement rate, however, is non-uniform across depth, and the recoverable resolution remains low near the native object plane; additionally this region exhibits strong artifacts. This effect is due to the depth-dependent sampling patterns and induced aliasing in light field imaging [19]. As the sample is naturally placed at the native object plane during the acquisition, i.e. in focus, the aliasing artifacts constitute a rather prominent problem in light field microscopy as expressed in [3,5,23], undermining the potential of the modality.
The sampling patterns and angular aliasing have previously been studied for light field imaging systems, like camera arrays [24][25][26][27]. However, there are fundamental differences between the aliasing in camera arrays and plenoptic devices [19], which must be acknowledged in order to address the cause of the artifacts in LFM deconvolution. Plenoptic devices avoid angular aliasing while introducing considerable spatial aliasing, since neighboring emitters in the scene are projected to pixels far apart on the sensor. Ng et al. [8] analyzed aliasing in refocused light fields, and Georgiev et al. [28] discussed the impact of the micro-lens array to sensor distance on the sampling rate in plenoptic cameras. In [19,29], the authors studied the depth-dependent sampling requirements in light field cameras.
In this work we study the depth-dependent sampling patterns of a light field microscope and analyze how they introduce aliasing, in order to understand the cause of the artifacts; to the extent of our knowledge, such an analysis has not been performed for the LFM before. We then derive depth-dependent anti-aliasing filters and propose a novel and efficient aliasing-aware deconvolution method for artifact-free 3D reconstruction.
Our work and [19] share the anti-aliasing priors idea in an interesting way. In [19], the authors use light field projections of the filtering kernels directly on the micro-images to ensure correct, non-aliased disparity maps. They then incorporate the estimated disparity in the light propagation model and proceed to recover the 2D radiance (all-in-focus image) from a light field camera image employing a variational Bayesian framework. This can be interpreted as an implicit accounting for the aliasing through the disparity prior. In our work, on the contrary, we derive the anti-aliasing filter kernels in the object space and explicitly apply them to the light field as a correction step of our iterative aliasing-aware deconvolution, which employs a smoothing expectation maximization scheme.
The analysis and the deconvolution scheme we propose apply to arbitrary plenoptic configurations. Hence, we also derive a generalized wave-based forward light propagation model able to characterize both original [1,22] and "focused" LFM setups. In previous works, the "focused plenoptic" camera [30,31] design was proposed to enhance the spatial resolution of the captured light field, compared to the original plenoptic camera [9], by manipulating the placement of the sensor with respect to the micro-lens array (MLA) such that the micro-lenses are focused on the objective lens. When coming to the light field microscope, due to the presence of the tube lens, manipulating the distance between the MLA and the camera sensor immediately affects the distance between the tube lens and the MLA such that the conjugate image of the native object plane may be in front or behind the MLA, creating a defocused field incident on the MLA, as opposed to the original LFM case, where the native image pane coincides with the MLA plane. Thus, although conflicting with the established "focused plenoptic" term, we find the term "defocused LFM" to better reflect generic LFM designs, as the recoverable resolution at a certain depth in the object space strongly depends on the extent of the defocus generated at the MLA plane.
We discuss the depth-dependent trade-offs in terms of recoverable lateral resolution when comparing various LFM configurations, and evaluate the proposed reconstruction algorithm on real experimental data to demonstrate high quality reconstructions, superior to previous work. As a teaser, Fig. 1(c) shows the reconstruction of the USAF 1951 resolution target to compare our proposed method (right) with the baseline method in [22] (center).
In summary, the contributions of the present work include: Analysis of the depth-dependent sampling patterns and aliasing in light field microscopy (section 2); a generalized wave-based o z = z) in front of the microscope objective has a conjugate image by the tube lens at z . Finally, the micro-lenses create micro-images at z , and the light reaches the camera sensor, producing a raw light field image. (b) Depth-dependent aliasing in LFM: The source points in the red group at depth z 0 in front of the microscope have completely overlapped images at the sensor plane. The points in the blue group, while being sampled at the same rate as the points in the red group, show partially non-overlapping images on the sensor as they are placed at depth z 1 . The points in the green group, on the other side, are also placed at z 1 ; however they are sampled at a higher rate and their images are fully non-overlapping. (c) Reconstruction of the USAF 1951 target: Left: a light field image of the USAF 1951 resolution target, acquired with our experimental LFM. Center: reconstructed target using the method in [22]. Specific aliasing artifacts are present. Right: artifact-free reconstruction using our aliasing-aware deconvolution method. light field point spread function to characterize defocused LFM designs (section 3); a novel and efficient aliasing-aware deconvolution method for artifact-free 3D reconstruction employing depth-dependent anti-aliasing filters (section 4); discussion and comparison of various LFM designs with respect to the recoverable resolution in the light field data and evaluation of the proposed deconvolution scheme on in-vivo biological samples and phantoms (section 5). Table 1 contains a summary of the symbols used in this paper together with their definitions.

The light field microscope
A light field microscope (LFM) is built by placing a micro-lens array (MLA) into the optical path of a conventional wide-field microscope [1,3]. Fig. 1(a) shows a ray diagram as an intuitive overview of the light propagation through a LFM. A source point at a depth z in front of the microscope objective has a conjugate image by the tube lens at z . The objective lens creates a virtual image of the object at z , which is not drawn here for the sake of clarity. We choose to represent z as f ob j + ∆z, since an object at depth f ob j (the focal length of the objective) is usually in focus in the wide-field microscope. In order to be consistent with the literature [3,22], we will call this depth, z = f ob j , the native object plane (NOP) or the zero plane of the LFM. Then ∆z represents an offset from the native object plane, and we will refer to this quantity when talking about depth in the subsequent sections. Finally, the micro-lenses create micro-images at z , and the light reaches the camera sensor, producing a raw light field image. Without loss of generality, Fig. 1(a) depicts a configuration where the conjugate image is formed in front of the MLA. However, throughout the paper, our derivations are valid for arbitrary configurations, i.e. they do not discriminate between defocused [30] and original plenoptic [1,8] light field imaging designs.

Light field microscope
A light field microscope (LFM) is built by placing a micro-lens array (MLA) into the optical path of a conventional wide-field microscope [1,3]. Figure 1(a) shows a ray diagram as an intuitive overview of the light propagation through a LFM. A source point at a depth z in front of the microscope objective has a conjugate image by the tube lens at z . The objective lens creates a virtual image of the object at z , which is not drawn here for the sake of clarity. We choose to represent z as f obj + ∆z, since an object at depth f obj (the focal length of the objective) is usually in focus in the wide-field microscope. In order to be consistent with the literature [3,22], we will call this depth, z = f obj , the native object plane (NOP) or the zero plane of the LFM. Then ∆z represents an offset from the native object plane, and we will refer to this quantity when talking about depth in the subsequent sections. Finally, the micro-lenses create micro-images at z , and the light reaches the camera sensor, producing a raw light field image.
Without loss of generality, Fig. 1(a) depicts a configuration where the conjugate image is formed in front of the MLA. However, throughout the paper, our derivations are valid for arbitrary configurations, i.e. they do not discriminate between defocused [30] and original plenoptic [1,8] light field imaging designs.

Depth-dependent sampling patterns of the LFM
We now proceed at investigating the sampling requirements of the LFM and deriving the depth-dependent quantities relevant for our algorithm, similar to the analysis of sampling patterns in plenoptic cameras [19]. Figure 1(b) is meant to build an intuition on the depth-dependence of the sampling patterns of the plenoptic microscope. The source points at z 0 in front of the microscope (circled in red) have completely overlapped images at the sensor plane. On the other side, the source points in the blue group have partially non-overlapping images on the sensor, while being sampled at the same rate as the points in the red group, but they originate from a different depth z 1 . The points in the green group, also at z 1 , are being sampled at a higher rate such that their images on the sensor are fully non-overlapping.
In order to characterize the depth-dependent nature of the sampling in light field microscopy, let us assume for now that the micro-lenses have very small apertures and behave like pinholes. Then we can approximate the MLA by an array of pinholes with spacing p ml . The in-camera light field at the MLA (pinholes in this context) should be band-limited with a bandwidth of f 0 = 1 2p ml in order to satisfy the Nyquist criteria [32]. Higher frequencies, outside this bandwidth, would be under-sampled by the pinhole array and appear aliased.
Since the sensor elements have a finite extent, we must look into what area the pixels effectively integrate over. Figure 2(a) illustrates how the image at the MLA scales to the actual image that forms under a micro-lens. For a clear visualization, we omit here the first part of the image formation and assume we have an image of an object at z in front of the objective formed at z by the objective lens. The tube lens further creates a scaled image at the conjugate image plane (dark blue), z . The image at the MLA (light blue) follows from tracking the chief rays. Finally, we pick a micro-lens and derive the micro-image behind it. By means of similar triangles, the size of the image under a micro-lens is the size of the image at the MLA, scaled by a factor: Here, d sens mla is the distance between the MLA plane and the sensor plane, and d mla tl is the distance between the tube lens and the MLA plane. The scaling amount γ z is depth-dependent, which means the actual area of the light field the sensor pixels integrate over varies with the object depth.
An interesting observation follows for telecentric microscopes (4f-systems), as in [3,22]: for these systems, although the magnification stays constant with object depth, the blur radius at the MLA (depicted in blue in Fig. 2(b) varies in extent with depth, Here, r tl is the effective tube lens radius. Consequently, the depth-dependent scaling factor γ z still applies, as we can write Please note the relation to the magnification factor in [19], or the amount of refocusing in [9].
If we now drop the pinhole array approximation and consider a micro-lens with finite aperture p ml , we have to take into account the additional blur they introduce. The depth-dependent blur  Fig. 3(a)), the NOP is sampled at the coarsest rate by the LFM which implies our resampling anti-aliasing filter has the largest radius at this object depth. As we move away from the zero plane, the LFM sampling rate increases and the anti-aliasing requirements become milder. under each micro-lens, depicted in red in Fig. 2 where r ml = p ml 2 is the radius of the micro-lens. We have now derived all the ingredients we need to characterize the non-aliasing requirements of the LFM.

Anti-aliasing filters
The band-limited assumption we made in the previous section for the pinhole approximation of the MLA means the acquired light field is the conjugate light field at z , convolved by a low-pass ideal (sinc) filter with cutoff frequency 1 2p ml . We define the sinc kernel radius as the first zero crossing of the filter, p ml . Then, as every micro-image is the projection of the conjugate light field onto the sensor, we can project the filter in the same way. Employing Eq. (1), the scaled filter kernel has a radius of γ z p ml .
When we take into account the finite micro-lens apertures p ml , the pixels effectively integrate over a larger area and the aliasing is reduced with the micro-lens blur b z . Then the filter size at the sensor, accounting for the micro-lens blur, is: In the case where the conjugate image forms at the MLA (z → d mla tl ), we have z → 0 and b z → ∞. However, the micro-lens blur actually converges to the size of the micro-image and thus we restrict the maximum filter radius to r ml : We now backproject the filter into the object space. For this we introduce the super-resolution factor, s ∈ Z, as defined in [22]. If we sample the volume at a rate of s times the lenslet resolution p ml , then the voxels are spaced by p ml Ms , where M is the objective magnification factor. Then the radius of our ideal filter kernel in pixels in object space is: Figure 2(c) illustrates the scaled pinhole filter radius together with the micro-lens blur radius and the final compensated anti-aliasing filter radius in pixels over an axial range [−100, 100]µm for an original plenoptic LFM configuration as in Fig. 3(a). An important observation here is that, as we move away from the zero plane, the LFM samples at a higher rate, imposing milder anti-aliasing requirements. Finally, we define h f w,z as the anti-aliasing normalized resampling filter. h f w,z is a depthdependent ideal low-pass filter and its kernel size at each depth is given by w obj z .
Since the ideal filter is of unit value for all the frequencies inside the band-limit, and zero outside, it has infinite extent. In practice, we need to use an approximate non-ideal filter kernel, aiming at optimizing the unity gain in the pass-band and zero gain in the stop-band. While there are extensive filter design choices [32], for all the experiments in this paper, we obtained great results using a Lanczos2 windowed version of the sinc kernel.

Generalized light field point spread function
In this section we propose a generalized forward light propagation model describing the optical system's impulse response for arbitrary LFM configurations (i.e. the light field point spread function, LFPSF).   In order to evaluate the wavefront incident on the MLA produced by a source point in front of the microscope, we employ Abbe imaging theory for 4-f optical systems [33]. We proceed to find the "focused" configuration for our scenario. In Fig. 4(b) FOP represents the focused object plane, which is the depth in the object space that is imaged exactly at the MLA by the tube lens. This plane is then located at offset ∆ NOP from the native object plane. If we introduce tl − f tl to be the signed distance between the MLA plane and the tube lens, we can write: where 1 M 2 is the axial magnification factor [33]. We now introduce the axial coordinate dof mla = f obj + ∆ NOP as the object space depth that is  from the FOP. Then we observe a defocused wavefront at the MLA plane: which is the Debye integral for circular lens apertures. deM = dof mla d mla tl is the demagnification factor, α ≈ arcsin(NA/n) is the maximum entrance angle of the objective aperture, n is the refractive index, λ is the wavelength of the monochromatic light we assume, P(θ) is the apodization function of the microscope, J 0 the zeroth order Bessel function of the first kind, and v, u are the normalized radial and axial optical coordinates respectively, given by: ∆z FOP represents the depth offset from the FOP for a source point, In order to stay consistent with the convention and for the clarity of the subsequent discussion, we will still refer to ∆z (o z = ∆z + f obj ) as the axial range of an object, via the following convenient substitution: ∆z FOP = ∆z + ∆ NOP (11) An immediate observation follows when ∆ NOP = 0, then dof mla = f obj and d mla tl = f tl . This is the original LFM configuration, and Eq. (9) is equivalent to the defocused PSF at the native image plane proposed in Broxton et al. [22].
Similar to us, in [2] the authors compute the wavefront at the MLA in a defocused LFM setup. They first model the PSF at the NIP as in [3,22] and then propagate the wavefront for a distance of ∆ MLA (the signed distance between the tube lens and the MLA plane) via Fresnel diffraction integral. In theory this is equivalent to our approach. In practice, however, FFT-based Fresnel propagation is implemented via its transfer function, which is a chirp function, and special sampling regimes have to be considered. This implies not only a computational overkill, but such requirements depend on the propagation distance and under-/over-sampling of the angular spectra introduce artifacts in the observation plane [34].

MLA transmittance
Having computed the light field at the plane immediately before the MLA, we account now for the effect of the MLA. The field U mla+ immediately after the MLA is given by: where the MLA transmittance function T is modeled by replicating the single lenslet transmittance in a tiled fashion as in [22]: with rep d,d the 2D replication operator and p ml the spacing between micro-lenses. t(x l , y l ) is the complex transmittance function of a lenslet with local lenslet coordinates, (x l , y l ): The exponential term is responsible for the phase change in the incident light, while P(x, y) represents the pupil function, where P(x, y) = circ p ml (x, y) for circular aperture lenslets or P(x, y) = rect p ml (x, y) for squared shaped lenslets. k = 2π λ is the wavenumber.

MLA to sensor light field propagation
We now address the propagation of the field from the MLA plane to the camera plane. Since we aim to model arbitrary distances between the MLA and the sensor, without restricting the d sens mla distance to satisfy the Fresnel number (paraxial assumption) [22], we use the more accurate Rayleigh-Sommerfeld diffraction solution [35] to predict the light field at the sensor plane: where (x s , y s ) are the image plane coordinates, F denotes the Fourier transform, and (f X , f Y ) are the spatial frequencies at the sensor plane. H rs is the Rayleigh-Sommerfeld transfer function:

Sampling requirements of the Rayleigh-Sommerfeld transfer function
The transfer function in Eq. (16) is a "chirp" function; it contains a phase function whose absolute value increases with the square of the frequency. Sampling such a function in practice is problematic since it is not bandlimited. In [34], the authors analyze the sampling requirements for the Rayleigh-Sommerfeld FFT-based propagation simulation and propose the following ideal sampling criteria: where L is the side length of the source plane and ∆x is the sampling interval such that f X ∈ [−1/2∆x, 1/2∆x] and ∆f X = 1/L. When this criteria is not met, and the transfer function in Eq. (16) is over-/under-sampled, the image might exhibit repeating patterns, aliasing and spike-like artifacts [34].
Since the source field U mla+ needs to match the sampling rate of the kernel H rs when implementing FFT-based propagation, and since the sampling requirements in Eq. (17) vary with the d sens mla distance, in practice, we implement a resampling strategy to satisfy these criteria.

F-number matching condition for defocused LFM setups
In order to ensure the micro-images optimally fill the sensor plane without overlapping when acquiring light field images, the effective image-side NA of the tube lens needs to match the effective NA of the micro-lenses. As depicted in Fig. 4(d), it is important to notice that a point source o(o x , o y , o z = dof mla ), generates the largest response (blur) behind a micro-lens. Conversely, as we move away from dof mla , the micro-lens blur radius, b z decreases. Thus, we only need to constrain the maximum blur radius, b dof mla to match the micro-lens radius r ml in order to ensure optimal non-overlapping sensor plane coverage. From Fig. 4(d) it quickly follows: where r tl represents the effective tube radius; the radius of the field distribution incident on the tube lens by a source point at dof mla in front of the microscope. In practice, we compute the r tl following the marginal rays: where z is obtained using the thin lens equation and d tl obj = f obj + f tl for 4f microscopes.
An immediate observation follows that when d mla tl = f tl and d sens mla = f ml , we have r tl = r obj (radius of the objective lens), and Eq. (18) is equivalent to M 2NA obj = f ml p ml , where M is the objective magnification and NA obj the objective numerical aperture. This is the f-number matching condition for the original LFM [1,22].
While violations of Eq. (18) result in either suboptimal sensor plane coverage or overlapping micro-images (see Fig. 4(e)), the LFPSF we derived in the current section allows for arbitrary d mla tl , d sens mla combinations and is consequently not limited to f-number matching configurations.

Aliasing aware 3D deconvolution
Having discussed the non-aliasing sampling requirements of the LFM and derived the generalized LFPSF, we now turn our attention to incorporating this prior knowledge into the reconstruction process of computing a 3D volume from a light field image.

Discretized imaging model
Given the raw noisy light field sensor measurements acquired m = (m j ) j ∈J by pixels j ∈ J (|J| = m) we seek to recover the fluorescence intensity at each discrete point in the volume which produced these measurements.
We represent the discretized volume v by a coefficient vector (v i ) i∈I with |I| = n. Note that the sampling rate in v is dictated by the super-resolution factor s defined in the previous section. We now denote the detection probabilities a ji = P(photon counted at sensor element j| emission occurred in voxel i) (20) Due to the low photon counts in fluorescence microscopy, we define the number of photons emitted at voxel i and detected by sensor element j as random variables z ji with z ji ∼ Poisson(v i a ji ), which we combine into the the iid random vector z = (z ji ) j ∈J,i∈I . Our measurements m = (m j ) j ∈J arise from z ji as m j = i∈I v i a ji , yielding the stochastic imaging model m ∼ Poisson(Av), where m denotes the light field measurement, v denotes the discretized volume we seek to reconstruct, and the operator A = (a ji ) j ∈J,i∈I describes the light field forward model, which is effectively determined by the discretized version of the LFPSF in Eq. (15). For each point in a fluorescent object the image intensity is given by the modulus squared of its amplitude [33]: where o(i) is the object space coordinate of voxel i and x s (j) is the coordinate of sensor pixel j.

Estimate-maximize-smooth algorithm
We now consider the estimation of v by maximizing the Poisson log-likelihood If we look at z as the complete version of the incomplete data m, the expectation maximization approach provides an iterative two-step scheme for increasing the likelihood of the current estimate v. In the first step, z is estimated by computing the conditional expectation E(z ij |m, v), and in the second step, the maximum likelihood estimate of v is found, starting from an initial guess v 0 :ẑ where q is the iteration count. This is the well known MLEM algorithm, and Eq. (25) also corresponds to the popular Richardson-Lucy [36] iterative update, which in matrix-vector notation reads: We now propose an additional straightforward step in which we filter the result of Eq. (24) and Eq. (25) using the depth-dependent anti-aliasing filters h f w,z that we derived in the previous section. Then the aliasing aware update scheme reads: where * represents the convolution operator. This smoothing step is computationally insignificant compared to the estimation and maximization steps.
The reconstructed v has a uniform lateral resolution across depths as imposed by the depth invariant discretization we choose; see the super-resolution factor s in the previous section. However, the non-aliasing sampling requirements of the LFM vary across depth, and the actual details that can be recovered depend on these patterns, among other factors.
Moreover, as our model does not incorporate explicit depth priors, information from one depth appears aliased when wrongly projected to another depth. This behavior is present from the first iteration of the Richardson-Lucy scheme, resulting in strong artifacts at the highly under-sampled depths where the process fails to converge. Thus, the resampling correction (by depth-dependent filtering) we propose is absolutely necessary.
The filtering in Eq. (27) can be interpreted as projecting v q+1 to the set of true solutions, which consists of frequencies below the bandwidth dictated by the LFM sampling requirements at each reconstructed depth.

Convergence of the proposed scheme
In order to show convergence of our proposed algorithm, we use the results of [37], in which a similar EMS algorithm with a smoothing kernel S is investigated. The authors in [37] demonstrate a modified (weighted) EMS algorithm with desirable convergence properties using a weighted smoothing kernel According to Lemma in Section 5.3 in [37], S and T will have approximately the same effect if S and W satisfy the three requirements: In our context, s i is the size of the voxel i, e i = j ∈J a ij and θ . For our smoothing kernel, h f w,z , the first two requirements are trivially fulfilled. Regarding the third requirement, as we use a uniform sampling of our volume, we have s i = s j ∀i, j. Also ∀i, j where h f w,z is non-zero, we have θ q i θ q j , as the effective sampling rate is higher than the kernel radius. The same argument holds for e i and e j , since the columns a i and a j act on regions of v at most as large as the effective sampling rate of the LFM. Thus we can conclude w i w j when S ji 0 and δ small.

Experimental results
All experiments in this work were performed with a custom-built LFM setup, configured as a 4-f system, combining a 0.5 NA with 20× magnification objective lens, and a tube lens with focal length f tl = 165mm. We used a f-number matching square-shaped aperture MLA with 150µm micro-lens pitch and 3000µm focal lengths. The pixel pitch of the sCMOS camera is 6.5µm, yielding a total of 23 × 23 pixels behind a micro-lens.
All the results we discuss in this section were reconstructed at sensor resolution, i.e. at a super-resolution factor s = 23, which translates to a uniform lateral resolution of 0.33µm. Note, this refers to the sampling rate we chose for rendering the volumes and has nothing to do with the actual details that can be recovered, which is the effective resolution of the LFM. We refer the interested reader to existing discussions on the subject [1,31,38].

Artifact-free deconvolution
In order to show the full potential of our method and to be able to use the method [22] as a baseline for comparison, the experiments in this section were done with the LFM in the original plenoptic configuration, i.e. d mla tl = f tl and d sens mla = f ml . Then the zero plane (∆z = 0µm) has a conjugate image exactly at the MLA and is the most under-sampled, exhibiting the most artifacts.
As the first experiment, a USAF 1951 resolution target was imaged at ∆z = 0µm using the LFM, see Fig. 1(c) (left) for the raw light field image. Figure 1(c) shows the reconstruction using the baseline method from [22] (center), and it is obviously riddled with the typical zero plane artifacts. Conversely, Fig. 1(c) (right) shows the reconstruction of the same light field image with our proposed aliasing-aware algorithm, exhibiting a natural appearance with no artifacts.
To further characterize the recoverable lateral resolution of our light field microscope configuration, we performed an analysis similar to [22]. We imaged and subsequently reconstructed the USAF target over a depth range of [−80, 80]µm, with 1µm steps. We then computed the contrast measure: C = I max − I min I max + I min (28) for each region of interest on the USAF target, from element 6.1, representing a spatial frequency of 64[lp/mm] to element 7.6 with spatial frequency 228[lp/mm]. In Fig. 5, we show, for both the baseline method [22] in (a) and our proposed method in (b), the lateral MTF by plotting the contrast of different line pair regions in the USAF reconstruction, over the axial position, together with example single plane reconstructed images. In (a), due to the presence of artifacts, the drop in resolution near the zero plane is not discernible, as artifacts are perceived as high contrast features. Conversely, in (b), the contrast plots match the expected resolution profile.
While in the state-of-the-art deconvolution [22], the aliasing artifacts make it difficult to quantify the lateral recoverable resolution in the highly affected axial ranges, when using our proposed EMS scheme, due to the smoothing step, the resolution can be reliably measured. However, although visually improved in quality, in terms of actual details, the effective resolution is ultimately limited by the sampling patterns and the PSF of the system. We encourage the reader to zoom into the example reconstructed USAF images in Fig. 5(a) and 5(b) for comparison.
As discussed in Section 3, when moving away from ∆z = 0µm, we need milder anti-aliasing filters to remove the artifacts while keeping the details in the reconstruction. This effect is illustrated in Fig. 6, which shows the reconstruction of an eyeball of a zebrafish larvae (5 days   [22] shows strong specific aliasing artifacts (red arrows) at depth planes close to the zero plane, while they fade out as we move further away from this plane. In comparison, our aliasing-aware deconvolution scheme completely removes all the artifacts. Figures 6(a) and 6(c) show the reconstruction with the baseline method in [22] and the artifacts are strongly present at depths close to the zero plane, while they fade out as we move further away from this plane, see for example the slice at ∆z = 25µm. In Figs. 6(b) and 6(d) we show the reconstruction of the same light field data using our proposed method. The depth-dependent filter radius is shown in Fig. 2(c); note here how the kernel radius drops as the artifacts fade away. Our deconvolution produces superior artifact-free results compared to the reference method without over-smoothing, as the depth-dependent filter is dictated by the sampling requirements of the LFM. Both reconstructions (Fig. 6) were obtained after 8 iterations of the corresponding update schemes. Figure 7 shows a reconstructed cardiomyocyte organoid labelled with the calcium dye Fluo4-AM. Such organoids are an emerging platform for clinical trials, enabling high-throughput studies (both logistically and ethically since they are not organisms), along with possible direct clinical applicability, since they are derived from human stems cells. Contrary to traditional cell-cultures they are extended across all three dimensions and they have interesting temporal dynamics like cell-signaling and, specifically for heart-organoids, movements.  [22] shows strong specific aliasing artifacts (red arrows) at depth planes close to the zero plane, while as we move away from this plane, the artifacts are less visible. Our aliasing-aware deconvolution method shows superior artifact-free results.
For the reference method [22], reconstruction artifacts are again strongly present (Figs. 7(a) and 7(c)) around ∆z = 0µm, while gradually fading out further away; see the slice at ∆z = 50µm. In the corrupted regions, a subsequent data analysis is not only troublesome, but rather unreliable. In comparison, our proposed method, specifically treats these artifacts via our depth-dependent resampling strategy (Figs. 7(b) and 7(d)). Again 8 iterations were performed for both methods.

Defocused LFM design
In this section we evaluate the defocused LFM setup. For this purpose, we place the micro-lens array at a distance d mla tl f tl from the tube lens, while keeping the tele-centricity as before (see

Figs. 3(b) and 3(c)
). The d sens mla distance will then follow from Eq. (18) to ensure the micro-images optimally cover the sensor plane. Figure 8 shows the reconstruction of a zebrafish eye (a different sample from the one in Fig. 6), over an axial range ∆z = [− 40,40]µm, from LF images acquired when d mla tl >f tl and d mla tl <f tl . In order to perform deconvolution on these images, we used the LFPSF we derived in Section 3. to describe the forward imaging model. In Fig. 8(b) we show the reconstruction using the classic Richardson-Lucy scheme, as in Eq. (26). While in the original LFM configuration the reconstructed samples show prominent aliasing artifacts around ∆z = 0µm due to the coarse sampling in this region, here, in the defocused LFM design we observe the artifacts pushed to the edges of our axial range. This effect is due to the displacement of the MLA from the native image plane, ∆ MLA , which effectively translates into a proportional axial shift of the depth dependent lateral sampling rates in the object space by ∆ NOP , see Eq. (8). This means an axial range ∆z in front of a defocused LFM setup is sampled in the same way the ∆z FOP range would be sampled by the original LFM setup with the same settings.  In Fig. 9, we imaged 1µm fluorescent beads in agarose with the LFM in the three configurations depicted in Fig. 3. Fig. 9(a) (top) shows the acquired LF image for the original plenoptic design (d mla tl = f tl ), together with the 3D reconstruction using our proposed method in Fig. 9(b) (top). Fig. 9(a) (middle) and (bottom) illustrate the acquired LF image for the defocused design with d mla tl > f tl and d mla tl < f tl , respectively, alongside the 3D reconstructions in Fig. 9(b) (middle) and (bottom). The red rectangle highlights two micro-spheres at the zero plane of the LFM (∆z = 0). While, in the original LFM setup (top) they are only reconstructed at the lenslet resolution, in Fig. 9 (middle) and (bottom) they appear better resolved. On the other side, the dashed arrow, for example, points to a sphere placed at the FOP plane in the defocused d mla tl > f tl (middle) case. While in this case it can only be recovered at lenslet resolution, in Fig. 9 (top) and (bottom) we observe it at higher resolution. Analogous discussion applies to the other beads. Just as we have seen for the fish eye in Fig. 8, while one LFM configuration performs well at spatially resolving certain depths in the axial range, it does so at the cost of other depths, which is also the case when imaging away from the zero plane with the original LFM. In order to extend the resolvable range in the 3-D reconstructions, such LFM configurations can be complementary, which supports and motivates the work towards multi-focus [39,40] or dual-camera [5,41] plenoptic setups.

Discussion
While we showed good results using a Lanczos2 windowed ideal filter (with depth-dependent cut-off frequencies) as the anti-aliasing resampling strategy, more advanced sampling-specific interpolation schemes should be considered to cope with the highly irregular sampling patterns of the LFM across depth, i.e. depth-dependent filter shapes in addition to the proposed depth-dependent filter sizes.
Although we improve the visual appearance of the 3D reconstructed objects by addressing the sampling artifacts, the effective lateral resolution is still limited by the depth-dependent sampling rate. We showed that defocused LFM designs better resolve the near native object plane range, Table 2 contains the system parameters for all the data sets we discuss in this paper, together with the relevant reconstruction settings. The Fish eye(>) and Fish eye(<) entries correspond to the d mla tl >f tl and d mla tl >f tl configurations in Fig. 8, respectively. Reconstructing the [−40, 40]µm axial range in both situations is equivalent, in terms of recoverable resolution, to reconstructing the [−80, 0]µm and [5,85]µm in the original LFM setup (see ∆z FOP column); effectively shifting the zero plane by ∆ NOP . This explains the strongest artifacts in Fig. 8 Fig. 8(c) we show the artifact-free deconvolution obtained using our Estimate-Maximize-Smooth scheme, and in Fig. 8(d) we illustrate z-slices of the reconstructed volumes every 5µm. We perceive better spatially resolved features in the range [0, 40]µm for the d mla tl >f tl setup, while in the d mla tl <f tl case this range appears blurred and the [−40, 0]µm range is resolved better.
In Fig. 9, we imaged 1µm fluorescent beads in agarose with the LFM in the three configurations depicted in Fig. 3. Figure 9(a) (top) shows the acquired LF image for the original plenoptic design (d mla tl = f tl ), together with the 3D reconstruction using our proposed method in Fig. 9(b) (top). Figure 9(a) (middle) and (bottom) illustrate the acquired LF image for the defocused design with d mla tl >f tl and d mla tl <f tl , respectively, alongside the 3D reconstructions in Fig. 9(b) (middle) and (bottom). The red rectangle highlights two micro-spheres at the zero plane of the LFM (∆z = 0). While, in the original LFM setup (top) they are only reconstructed at the lenslet resolution, in Fig. 9 (middle) and (bottom) they appear better resolved. On the other side, the dashed arrow,  Table 2) version of the original LFM; the zero plane behavior is now appearing at the FOP plane. (c) The artifact-free reconstruction using our aliasing-aware deconvolution method. (d) Lateral slices through the reconstructed volume. The two defocused LFM configurations demonstrate higher resolved features at complementary axial ranges; marked with smileys. for example, points to a sphere placed at the FOP plane in the defocused d mla tl >f tl (middle) case. While in this case it can only be recovered at lenslet resolution, in Fig. 9 (top) and (bottom) we observe it at higher resolution. Analogous discussion applies to the other beads. Just as we have seen for the fish eye in Fig. 8, while one LFM configuration performs well at spatially resolving certain depths in the axial range, it does so at the cost of other depths, which is also the case when imaging away from the zero plane with the original LFM. In order to extend the resolvable range in the 3-D reconstructions, such LFM configurations can be complementary, which supports and motivates the work towards multi-focus [39,40] or dual-camera [5,41] plenoptic setups.

Discussion
While we showed good results using a Lanczos2 windowed ideal filter (with depth-dependent cut-off frequencies) as the anti-aliasing resampling strategy, more advanced sampling-specific interpolation schemes should be considered to cope with the highly irregular sampling patterns of the LFM across depth, i.e. depth-dependent filter shapes in addition to the proposed depth-dependent filter sizes.
Although we improve the visual appearance of the 3D reconstructed objects by addressing the sampling artifacts, the effective lateral resolution is still limited by the depth-dependent sampling rate. We showed that defocused LFM designs better resolve the near native object plane range, while sacrificing resolution at other depth ranges by axially shifting the sampling patterns; similar to imaging away from the zero plane with the original LFM design. In this regard, the LFM in any of the configurations, introduces depth-dependent trade-offs in terms of resolution by design.
In order to improve on these limitations, hybrid systems have been proposed, combining a light field microscope with a wide-field acquisition [42], or splitting the optical path and image with complementary focused LF systems [5,41]. Such dual-/multi-camera setups aim at achieving highly uniform sampling and subsequently isotropic resolution.
Alternatively, employing an array of lenslets with mixed focal lengths introduces more irregularity in the sampling patterns and allows increased depth of field [39,40]. [13] also demonstrates higher lateral resolutions over extended depth of field by introducing phase masks into the optical path of the LFM and thus creating highly variant PSFs for different depths. Interestingly, in [29], the authors discuss denser and more uniform sampling patterns in light field photography by incorporating lens aberrations and misalignment into the imaging model.
In order to improve the effective resolution, all of these approaches target, in one way or another, the depth-dependent nature of light field sampling. We believe, however, that a setup-specific anti-aliasing resampling strategy is still necessary in deconvolution schemes for such approaches.
Finally, now that we can eliminate the troublesome artifacts with our proposed method, it becomes feasible to improve the deconvolution method by including other image priors, such as edge-enhancing regularization; previously, such regularization only enhanced the artifacts.

Conclusion
In this work we address one of the current challenges in 3D reconstruction of light field microscopy data, the aliasing artifacts. We perform an analysis of the aliasing-free sampling requirements of the LFM to derive depth-dependent anti-aliasing filters. We also derive a generalized wave-based LFPSF to propose a novel aliasing-aware deconvolution scheme that applies to arbitrary LFM designs. We compare the capabilities of the original and defocused LFM designs in terms of recoverable lateral resolution at various axial ranges and demonstrate the superior quality reconstruction performance of our method using experimental data from phantoms and in-vivo biological samples.

Implementation and datasets
The example datasets shown in section 5. together with the implementation of the methods described in this paper are available as part of our 3D reconstruction framework for light field microscopy oLaF, available at: https://gitlab.lrz.de/IP/olaf.