## Abstract

Realizing both high temporal and spatial resolution across a large volume is a key challenge for 3D fluorescent imaging. Towards achieving this objective, we introduce an interferometric multifocus microscopy (iMFM) system, a combination of multifocus microscopy (MFM) with two opposing objective lenses. We show that the proposed iMFM is capable of simultaneously producing multiple focal plane interferometry that provides axial super-resolution and hence isotropic 3D resolution with a single exposure. We design and simulate the iMFM microscope by employing two special diffractive optical elements. The point spread function of this new iMFM microscope is simulated and the image formation model is given. For reconstruction, we use the Richardson-Lucy deconvolution algorithm with total variation regularization for 3D extended object recovery, and a maximum likelihood estimator (MLE) for single molecule tracking. A method for determining an initial axial position of the molecule is also proposed to improve the convergence of the MLE. We demonstrate both theoretically and numerically that isotropic 3D nanoscopic localization accuracy is achievable with an axial imaging range of 2um when tracking a fluorescent molecule in three dimensions and that the diffraction limited axial resolution can be improved by 3–4 times in the single shot wide-field 3D extended object recovery. We believe that iMFM will be a useful tool in 3D dynamic event imaging that requires both high temporal and spatial resolution.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Fluorescence optical microscopy is one of the most widely used imaging tools in molecular and cell biology because of its non-invasive and high biochemical labeling specificity for *in vivo* measurements. However, it has two fundamental limitations. First, it is too slow to capture many 3D dynamic events in a snapshot manner, due to the long acquisition time required to sequentially z scan the focal plane by moving either the object stage or objective lens, e.g. as in scanning confocal microscopy. Second, the axial spatial resolution (roughly 500–700 nm) is substantially lower than the lateral one (roughly 200–300 nm) because of the limited collecting angles of an objective lens. In a single objective configuration, the 3D intensity point spread function (PSF) features an elongated focal spot along the optical axis [1] and its optical transfer function (OTF) suffers from the well-known “missing cone” problem along the axial direction [1], posing a particular challenge in isotropic 3D microscopy imaging.

To overcome the 3D imaging speed limitation, multifocal plane microscopy (MUM) [2] was recently proposed to allow simultaneous acquisition of multiple focal planes in a single exposure time of a detector. Initially, this was implemented by using multiple beam splitters and detectors. In that case, each detector is placed at a specific distance from the tube lens and images a distinct focal plane within the sample [2, 3]. Currently, this setup is capable of imaging up to 4 distinct planes [4, 5] and becomes increasingly bulky if more planes are imaged since it requires one detector per focal plane. Another approach, called multifocus microscopy (MFM), utilizes a multi focus grating (MFG) in the Fourier plane to image 9 or 25 focal planes on a single detector [6–9]. In this method, a regular microscope is augmented by a 4f system, with a MFG placed in the Fourier plane. The MFG splits and diffracts the light into several beams traveling in different directions, and each beam is focus shifted so that different depth layers in a 3D volume can be projected onto different sub-regions of the imaging plane simultaneously at the expense of the lateral field [8]. To correct chromatic aberration (CA), a chromatic correction grating (CCG) and a multifacet prism are placed after the MFG [8]. Recently, we simplify the system by abandoning the CCG and prism. Instead, we place a narrow bandwidth filter in front of the detector to mitigate CA [10]. Since their advent, MUM and MFM have attracted significant interest in applications that involve the investigation of 3D dynamic samples. Examples include the tracking of a single molecule in three dimensions [3, 4, 7, 8, 11], or the detection of a dynamic process in the thick samples [5, 8, 12, 13].

However, like a conventional focal scanning microscope, both MUM and MFM also suffer from the anisotropic 3D resolution problem due to the limited collecting angles of the single objective lens. For instance, in MUM single-particle tracking, the axial localization accuracy is relatively worse than the lateral one [3, 4]. One scheme for achieving isotropic 3D super resolution with single fluorescent molecules is PSF engineering. For example, the PSF of the microscope has been engineered to have 2 rotating lobes where the angle of rotation depends on the axial position of the emitting molecule (so called double-helix PSF) [14]. However, this method is not demonstrated to achieve axial super resolution for 3D extended object though it works well for single molecule localization. Another common solution to achieving the axial super resolution and hence isotropic 3D resolution is to increase the solid aperture angle of the microscope by using two opposing lenses [1]. In the dual objective lens microscopy configuration, the fluorescent light emitted from the same molecule is collected through both lenses simultaneously, while interfering with each other on a common detector plane. The interference pattern contains high resolution information about the axial position of the molecule, which is otherwise truncated by a single objective lens resulting in poorer axial resolution in a single lens imaging. Typically, the axial resolution in the dual objective lens interferometry detection scheme attains 3 to 4-fold improvement over that obtained by a single lens. In addition, the axial resolution can be further improved by 1.5 to 2-fold if the excitation light is also split into two paths and coherently superposed to form axially varying standing wave (SW) illumination at the sample [15]. The coherent synthesis of two lenses illumination or/and detection is essentially the core concept in 4Pi point scanning microscopy, wide-field focal scanning image interference microscopy (I^{2}M) [1], and incoherent illumination interference image interference microscopy (I^{5}M) [1]. The dual objective interferometric detection scheme has also been adopted in localization-based super resolution imaging to achieve better axial localization precision. Examples include interferometric photoactivated localization microscopy (iPALM) [16], 4Pi single marker switching microscopy (4Pi-SMS) [17], and rapid interferometric particle tracking in 3D [18].

In this paper, we design and simulate a new microscope to increase 3D imaging speed and axial resolution simultaneously, for the first time. Specifically, we introduce an interferometric multifocus microscopy (iMFM) system, a combination of MFM with two opposing objective lenses that provides significantly improved axial and hence isotropic 3D resolution with a single shot. We design and simulate the proposed iMFM by employing two MFGs in the associated Fourier planes. We show that the use of two diffractive optical elements (DOEs) is necessary and key for iMFM to make sure the conjugate spherical wavefronts emitted from the same molecule are diffracted into the same diffraction order on the detector and hence interfere effectively after passing through two lenses. Both monochromatic and polychromatic PSFs of this new microscope are simulated and the image formation model is provided. We show by simulation that the proposed iMFM is capable of recording multiple interferometric focal planes simultaneously in a single shot, which contains axial super resolution information. We demonstrate the iMFM microscope with two simulated applications: single molecule tracking and 3D extended object recovery. The Fisher information matrix (FIM) and the Cramér-Rao lower bound (CRLB) of iMFM PSFs are theoretically calculated. A method to determine an initial axial position of the molecule is proposed to improve the convergence of the MLE localization algorithm in iMFM. We demonstrate both theoretically and numerically that isotropic 3D nanoscopic localization accuracy is achievable with an axial imaging range of 2um when tracking a fluorescent molecule in three dimensions and that the diffraction limited axial resolution can be improved by 3-fold in 3D extended object recovery with a single exposure in the proposed iMFM, significantly speeding up the acquisition process of conventional dual-objective scanning interferometric microscopes, i.e., I^{2}M [1].

The remainder of the paper is organized as follows. In section 2, we describe the theory and principle of the proposed iMFM system, with an emphasis on its PSF derivation and simulation. In sections 3 and 4, we describe two applications of iMFM: single molecule tracking and 3D extended object imaging, respectively. Both theoretical analysis and numerical reconstructions of the synthetic 3D images are given to demonstrate the superiority of the proposed iMFM with respect to axial resolution improvement. Conclusion is given in section 5.

## 2. iMFM principle

#### 2.1. Challenge of combining I^{2}M and MFM

A straightforward method to design an iMFM microscope might be to simply follow a dual-objective image interference microscopy (I^{2}M) with an MFM stage, as shown in Fig. 1(a). In such a design, a regular I^{2}M configuration is augmented by a 4f system (relay lens system L3 and L4) with a multi-focus grating (MFG) in the Fourier plane. The emitted light from the same fluorescent point simultaneously enters both upper and lower objectives with conjugate wavefronts, and then is coherently added to form the interference pattern on a common beam splitter (BS) plane, as in I^{2}M. However, when these two conjugate wavefronts emitted from the same fluorescent point propagate through the same MFG, they will be diffracted and focused into the positive and the negative diffraction orders on the detector [6, 8]. Looking at the color-coded focal planes shown in Fig. 1(a), we can see that the MFG will incorrectly interfere the orange focal plane from the top arm of the interferometer with the neon green focal plane from the bottom arm of the interferometer, and the yellow with the red, the green with the purple, and so on until the neon green focal plane of the top arm is interferometrically combined with the orange focal plane of the bottome arm. Therefore, in this naïve I^{2}M + MFM configuration, the MFG prohibits the correct interference to occur on the detector.

#### 2.2. Proposed iMFM framework

To overcome the above problem, our design of iMFM employs two MFGs in the respective Fourier planes, as shown in Fig. 1(b). This iMFM configuration features dual-objective detection scheme with an additional opposing objective lens. Each detection arm is a conventional MFM, which consists of an objective lens, a tube lens, and a 4f system (relay lens system L1 and L2) with an MFG in the Fourier plane. According to MFM principle [8], each objective lens images multiple focal planes into an array of differently focused tiles on the detector within a single exposure of time.

The key difference between iMFM [Fig. 1(b)] and the naïve I^{2}M + MFM [Fig. 1(a)] is that the emitted light is first split into multiple tiles and then coherently added on the common BS plane in iMFM. Hence, if two MFGs of opposite focal shifts are employed in the dual detection arms, the conjugate wavefronts from the same fluorescent point will be diffracted and focused on the same tile, and therefore interfere on the same area of the BS and the detector. That is, the complementary MFGs in each arm place like-colored focal planes of the sample onto the same tile on the detector. Therefore, the proposed iMFM is capable of producing a focal stack of interferometric 2D images simultaneously on a single detector with a single exposure without focal scanning. Note that even though we use two MFGs in the proposed iMFM setup, those two MFGs are actually identical. Because the MFG with an opposite focal shift can be simply obtained by rotating from one another by 180 degrees [see Fig. 2].

To correct CA, a CCG and a multifacet prism can be added at each detection arm of the proposed iMFM system shown in Fig. 1(b) by following the CA correction procedures outlined in [8]. In this paper, we consider both chromatic aberrated iMFM system and un-aberrated iMFM system. The reasons for proposing the chromatic aberrated iMFM system are twofold: (1) it allows us to simply the system complexity and cost by abandoning CCGs and prisms; (2) it allows to see the effect of the CA on the iMFM PSF and reconstruction quality.

#### 2.3. iMFM PSF and simulation

In this section, we will mathematically derive the iMFM PSF by assuming uniform excitation light at the sample. Our derivation is an extension of single-objective MFM [8] to the dual-objective iMFM. Suppose a single point source (0, 0, *z*) is sandwiched between two opposing objectives, where *z* is the displacement of the point from the sample focal plane. Upon fluorescence emission, the upper (obj1) and lower (obj2) objective lenses in Fig. 1(b) respectively produce complex wavefronts on the detector plane (*x, y*) as:

*x, y*) are spatial coordinates at the detector plane,

*λ*is the fluorescence emission wavelength, ℱ denotes the 2D Fourier transform,

*g*

_{1}and

*g*

_{2}denote MFGs placed in the Fourier planes of the two objective lenses, respectively, with spatial coordinates of (

*x*).

_{g}, y_{g}*f*

_{1}and

*f*

_{2}are complex wavefronts in the upper and lower objective Fourier plane caused by the point source (0, 0,

*z*) propagation, respectively. Let ℱ {

*f*

_{1}(

*x*; 0, 0,

_{g}, y_{g}*z*;

*λ*)} =

*p*

_{1}(

*x, y*; 0, 0,

*z*;

*λ*) and ℱ {

*f*

_{2}(

*x*; 0, 0,

_{g}, y_{g}*z*;

*λ*) =

*p*

_{2}(

*x, y*; 0, 0,

*z*;

*λ*), where

*p*

_{1}(

*x, y*; 0, 0,

*z*;

*λ*) is the coherent spread function (CSF) of a point source (0, 0,

*z*) in a single lens imaging under uniform illumination and has the following form in the focal region of a high numerical aperture (NA) objective of circular aperture under the scalar and Debye approximations [19]:

*A*is a constant,

*α*= sin

^{−1}(NA/

*n*

_{0}is the semi-aperture angle of the single objective lens, in which

*n*

_{0}is the index of refraction,

*k*= 2

*πn*

_{0}/

*λ*is the wave number,

*J*

_{0}is the zeroth order Bessel function of the first kind, and $\rho (x,y)=\sqrt{{x}^{2}+{y}^{2}}$ denotes the radial coordinate position on the detector plane. Note that $\sqrt{\mathrm{cos}\theta}$ is an apodization function for a high NA objective under the Abbe sine condition. The CSF modeling with a low NA objective can be derived from Eq. (3) by further assuming the paraxial approximation. The more details on such derivation and discussion can be found in [19]. In dual-objective detection, we have the following relationship [20, 21] because of opposite propagation directions of the emitted light into both upper and lower objectives.

For monochromatic light, the detection intensity PSF is the square of the coherent addition of spherical wavefronts of two opposing objective lenses, expressed as

*g*

_{1}and

*g*

_{2}are uniformly equal to one inside the circular pupil plane), Eq. (5) reduces to the detection PSF of I

^{2}M or I

^{5}M [1,20,21], which records interferometric information from one focal plane at a time, and then assembles a 3D interferometric stack from sequentially refocused 2D images.

In our iMFM configuration shown in Fig. 1(b), the entire 3D interferometric focal stack can be recorded in a single shot by placing two MFGs of opposite focal shifts in the Fourier planes of dual objectives. This is because the MFG can be designed to split the emitted beam into an array of (2*M* + 1) × (2*N* + 1) differently focused tiles on the detector at one exposure time [8,9], where *M* and *N* are non-negative integers, denoting the highest horizontal and vertical diffraction orders detected by the senor. The focal shift property of the MFG is achieved by imposing a geometrical distortion on a regular periodic grating pattern (with periods *d _{x}* and

*d*in the

_{y}*x*and

*y*directions, respectively) placed in the Fourier plane. This geometrical distortion in return introduces an order-dependent defocus phase shift in the detector plane. Mathematically, the MFG equation can be described as (see Appendix for more details):

*m, n*) such that $\sum}_{m=-M}^{M}{\displaystyle {\sum}_{m=-N}^{N}{w}_{m,n}^{2}\le 1$ is the total efficiency of MFG,

*z*= (

_{mn}*m*+

*Bn*)Δ

*z*is a focal shift at diffraction order (

*m, n*), in which Δ

*z*is a predefined focal step between two adjacent focal planes and

*B*= 2

*M*+ 1,

*λ*is the central emission wavelength used for MFG distortion calculation,

_{c}*δ*(

*x, y*) is a Dirac delta function, and (

*mx*

_{0},

*ny*

_{0}) = (

*mf λ*/

_{c}*d*/

_{x}, nf λ_{c}*d*) is the center position of diffractive order (

_{y}*m, n*) on the camera plane under the paraxial approximation for the emission wavelength

*λ*. Note that the paraxial approximation holds here because the focal length

_{c}*f*of the relay system lens L4 in Fig. 1(b) is much larger than the sensor size. The above equations indicate two distinct properties of the MFG: (1) light path splitting into an array of (2

*M*+ 1) × (2

*N*+ 1) diffraction orders (or tiles) indicated by Dirac delta function due to the periodic property of the MFG, and (2) order-dependent phase shift indicated by the exponential phase function due to the distortion of the MFG. Note that if a CCG and a prism are added in each arm of the iMFM system to correct CA, the modeling of the combined effect of the MFG and an optical correction module can be obtained by setting

*λ*/

*λ*= 1 in Eqs. (6) and (7), indicating that different wavelengths of the light focus at the same position without dispersion.

_{c}If the wavefront defocus phase from the out-of-focus plane, exp (*ikz* cos *θ*) as shown in Eq. (3), is compensated (or corrected) by the MFG order-dependent defocus phase, exp(−*ikz _{m}*,

*cos*

_{n}*θλ/λ*) as shown in Eq. (6), it is possible to simultaneously focus the light originating from an in-focus plane and multiple out-of-focus planes onto distinct lateral diffraction orders, and therefore form multi-focus images on a single 2D camera within one exposure time. This is the principle of the single lens MFM [8]. In our dual-objective iMFM configuration, the light originating from the same emitted point source has conjugate phase after passing through two opposing objectives. Therefore, in order to diffract these two conjugate spherical wavefronts to the same diffraction order on the detector, two MFGs

_{c}*g*

_{1}and

*g*

_{2}with opposite focal shifts are needed and have to be placed in the Fourier planes in the iMFM microscope, as shown in Fig. 1(b). The use of two phase masks to produce the correct interference between two coherent beams has also been recently reported in interferometric rotating double-helix (DH) PSF engineering [22].

According to the convolution theorem, by substituting Eqs. (3), (4), (6) and (7), Eqs. (1) and (2) become

*m*,

*n*). Note that in Eq. (10) a summation operator is taken outside a square operator because of the non-overlapping between each tile image by placing a physical field stop at the intermediate image plane in MFM [10].

The polychromatic intensity PSF is an integration of the monochromatic PSF over the emission spectrum Δ*λ*:

*z*,

_{m}*/*

_{n}λ*λ*for emission wavelength

_{c}*λ*, and more importantly, features axial interference pattern, which contains high frequencies and high resolution information. Note that by setting

*λ*/

*λ*= 1, Eq. (12) will become an analytical formulations of the PSF of the chromatic corrected iMFM system.

_{c}#### 2.4. Simulation parameters and results

To verify our iMFM system and make a practical comparison between MFM and iMFM PSFs, we have simulated the detection PSF for a numerical aperture NA = 1.27 water immersion objective, with refractive index of 1.338, the magnification $\widehat{M}=60$ and the emission central wavelength *λ _{c}* = 620nm. For simplicity, the focal lengths of all the relay lenses (L1-L4) are set to be 200mm, producing a unit magnification image relay. Note that in practice, the focal lengths can vary to introduce the non-unit magnification image relay if necessary. Based on these parameters, we designed a binary phase only (0 or

*π*) MFG1 [shown in Fig. 2(a)], which produces 3 × 3 focal shift images with a focal step of Δ

*z*= 250nm. The period of the MFG1 is

*d*=

_{x}*d*= 56

_{y}*µ*m with a pixel size of 1

*µ*m, and the diameter of MFG1 is set to be the same as that of the pupil aperture, defined by 2

*f*NA = 8467

_{obj}*µ*m. The unit cell pattern of the MFG1 is optimized for a maximum diffraction efficiency by using the iterative Gerchberg-Saxton (GS) algorithm [23]. As a result, the diffraction efficiency ${w}_{m,n}^{2}$ for each diffraction order is [7.13, 7.79, 7.01, 7.70, 8.57, 7.70, 7.01, 7.79, 7.13] percent with a total efficiency of 67.83%, which is close to the theoretical maximum efficiency [24]. Note that the higher uniformity of the intensity distribution between tiles can be obtained by using multiple random initial values for the GS algorithm and picking the optimal one. Also, the total efficiency can be improved by using a multi-phase MFG [24] or binary-phase MFG with phase different than

*π*[25]. The geometric distortion is then imposed on MFG1 to create a proper focal step of Δ

*z*= 250nm between consecutive diffraction orders for the emission central wavelength

*λ*. More details of the principle of MFG design can be found in [8,9]. The other MFG with opposite focal shift was generated by rotating MFG1 by 180 degrees [Fig. 2(b)]. The detection PSFs were numerically simulated in the size of 1500 × 1500 × 200 voxels with the voxel size of 80nm × 80nm × 10nm for both MFM and iMFM. The monochromatic unaberrated iMFM detection PSF is computed from Eq. (5) by setting

_{c}*λ*=

*λ*. For a polychromatic PSF in the presence of chromatic aberrations (CA), a 10nm emission filter is considered and the PSF is computed by integrating Eq. (5) over the emission bandwidth of 10nm. Note that all the parameters shown here are chosen to be similar to those reported for previous experimental MFM implementations [10].

_{c}The *xy, xz, yz* cuts and 1D axial profile of monochromatic and polychromatic MFM and iMFM PSFs are shown in Fig. 3. Note that for the axial PSF profile (right column of Fig. 3), the peak intensity of each tile is normalized relative to the tile’s intensity distribution ratio. Figure 3 verifies that both monochromatic and polychromatic iMFM PSFs are composed of 3 × 3 focal shift interferometric tile-PSFs, as indicated by Eqs. (10) and (12), respectively. When CA exists in the presence of 10nm bandwidth emission, each non-central tile-PSF has a dispersion along *x* or/and *y* directions, causing a huge loss of the peak intensity except for zero diffraction order, as shown in the 1D axial profile in Figs. 3(c) and 3(d). However, full width at half maximum (FWHM) along *z* direction remains almost the same, with about 190nm for iMFM PSF and 620nm for MFM PSF, with 3.3-fold decrease of the diffraction spot size along the axial direction.

#### 2.5. Image formation model and reconstruction in iMFM

The proposed dual-objective iMFM allows single-shot multifocal interferometry detection, mapping multiple interferometric focal planes of 3D sample volume *o*(*x, y, z*) simultaneously onto a 2D image plane *I*(*x, y*) in one exposure time without translating the sample. In this case, the recorded intensity is given by

*h*

_{iMFM}is the iMFM monochromatic or polychromatic PSF,

*b*is homogeneous background noise, and $\widehat{n}$ denotes additive Gaussian noise. For compaction, Eq. (13) can also be written in a matrix-vector form as where

**I**is an ℳ × 1 column vector in which ℳ = ℳ

*ℳ*

_{x}*is the number of pixels of the recorded image,*

_{y}**o**is an $\mathcal{N}\times 1$ column vector in which $\mathcal{N}={\mathcal{N}}_{x}{\mathcal{N}}_{y}{\mathcal{N}}_{z}$ is the number of voxels of the 3D unknown object to be recovered,

**A**is the $\mathcal{M}\times \mathcal{N}$ sensing matrix computed from iMFM’s normalized PSF

*h*

_{iMFM}, whose integral is equal to one. Additive Gaussian noise is ignored here because both MFM and iMFM are photon budget-limited detection schemes due to light splitting property of the MFG, and thus, the Poisson noise dominates [26].

In order to take the Poisson noise into account and remove the out-of-focus blur, we use the Richardson-Lucy (R-L) deconvolution algorithm [27] for iMFM 3D extended object recovery. To suppress noise amplification, total variation (TV) regularization is used. In addition, the non-negativity constraint is applied due to the non-negative nature of fluorescent light so as to restrict the set of possible solutions. These constraints are very helpful when we recover more 3D information from fewer 2D measurement data in MFM, i.e. when $\mathcal{M}<\mathcal{N}$ [10, 28]. In R-L deconvolution, the following optimization is performed [27]

*p*denotes the pixel coordinate in the captured image,

*λ*

_{TV}is the regularization parameter, and TV (

**o**) = Σ

*|∇*

_{s}**o**(

*s*)|, in which

*s*denotes the voxel coordinate in the

**o**. A solution to the optimization problem in Eq. (15) can be found by the following iteration [27]

*k*denotes the iteration number, the divisions are element wise,

**A**is the transpose of

^{t}**A**, and

*div*stands for the divergence operator. We notice that the denominator of Eq. (16) may become zero or negative due to a large value of

*λ*

_{TV}. To prevent this from happening, we set the negative values or “not a number” (NAN) values to zero at each iteration step. The algorithm terminates when the difference between two consecutive values of the cost function is smaller than a predefined threshold.

In single particle tracking, the 3D space consists of a single point with varying 3D positions over time. We model the single point as $o=\delta (\mathit{\theta}-\widehat{\mathit{\theta}})=\delta (x-{\widehat{x}}_{0},y-{\widehat{y}}_{0},z-{\widehat{z}}_{0})$ and use MLE to recover its 3D position $\widehat{\mathit{\theta}}=({\widehat{x}}_{0},{\widehat{y}}_{0},{\widehat{z}}_{0})$. In MLE, the following optimization is performed:

*f mincon*function.

We remark that the optimization problem in Eq. (17) for MLE localization is non-convex, and therefore it is very sensitive to the initial value. To improve the accuracy of MLE, we propose a method to determine the initial axial position of the single point that is imaged by the iMFM microscope next in section 3.

## 3. Isotropic 3D localization for iMFM single molecule tracking

#### 3.1. Fisher information matrix for the iMFM PSF

The Fisher information matrix (FIM) measures the sensitivity of an observation (e.g., iMFM PSF) to changes of the parameters to be estimated (e.g., 3D position of a single molecule). The model for calculating the FIM for iMFM is the same as that for MFM [3,26]. For each tile image, the photon detection is an independent Poisson process [26]. Therefore, the total FIM of an iMFM PSF is the sum of the FIM for each tile PSF, and can be expressed as a 3 × 3 matrix

*i, j*∈ [

*x, y, z*],

*N*is the number of pixels for each tile image,

_{p}*b*is the homogeneous background photons per pixel, ${\widehat{h}}_{m,n}(p)$ is a normalized tile-PSF at diffraction order (

*m, n*), and ${N}_{m,n}={w}_{m,n}^{2}2{N}_{\text{tot}}\gamma $ denotes the number of photons collected by the tile-PSF, where

*N*

_{tot}denotes the total number of photons collected by each objective lens and

*γ*= 0.5 denotes the photon loss ratio at the beam-splitter (BS). Then CRLB, which bounds the variance of the localization estimation

*σ*

^{2}, can be calculated by taking the diagonal elements of the inverse of the FIM as

It can be seen from Eq. (19) that the Fisher information or the localization precision can be improved by increasing the derivative of the PSF, i.e., $\partial {\widehat{h}}_{m,n}/\partial i$. In our dual-objective iMFM, each tile-PSF has a larger z-derivative due to its axial interference pattern [see Fig. 4(a)], and therefore greater axial differential information content, causing 3 to 4-fold improvement in axial localization [29] compared with single lens MFM. In addition, the simultaneous multifocal detection of iMFM leads to almost uniformly high combined differential information content along the large depth range [see Fig. 4(b)], while conventional dual-objective detection in a single channel suffers from non-uniform localization along *z* due to zero z-derivative at its PSF intensity nodes [29]. Figure 4 shows xz-cuts of the square of z-derivative of the normalized MFM and iMFM PSFs.

According to Eq. (19), the resolution can also be improved by increasing the number of the collected photons, i.e., *N _{m}*,

*. The dual-objective detection collects twice the number of photons compared with single lens imaging, and therefore improves the resolution by a factor of $\sqrt{2}$ in all three dimensions when the background*

_{n}*b*is small [30]. However, the light efficiency of our dual-objective iMFM configuration shown in Fig. 1(b) is similar to that of single objective MFM in that there is a factor of two light loss in the BS. To avoid the light loss, second detector can be introduced and placed at the output of the BS [29–31]. The axial resolution improvement is demonstrated in Fig. 5(left), in which we plot the theoretical localization precision

*σ*along 2

_{z}*µ*m axial range calculated from Eqs. (19) and (20). In the calculation,

*N*

_{tot}= 2500 detected signal photons per objective lens and

*b*= 10 background photons per pixel are considered, which are typical values observed in single molecule experiments [32–34] and remain the same throughout the paper unless stated otherwise. The parameter values of the microscope and MFG are set to be the same as those when we simulated the PSF, as described in section 2.3. The results suggest that the monochromatic iMFM PSF provides an average theoretical localization precision of (

*σ*) = (16.6nm, 16.6nm, 11.2nm) [red line] over the imaging range of 2

_{x}, σ_{y}, σ_{z}*µ*m for 2500 total signal photons per objective and 10 background photons per pixel. Compared to the single lens monochromatic MFM PSF with localization precision of (18.0nm, 18.0nm, 44.8nm) [blue line], the lateral localization precision $\sqrt{{\sigma}_{x}{\sigma}_{y}}$ is improved by a factor of 0.08 (this small improvement probably results from the redistribution of the light due to the interference), and the axial one by a factor of 4. We observe an oscillation of the z-localization precision along

*z*direction for iMFM. This is an inherent problem for the dual-objective interferometric microscope due to the zero z-derivative at its PSF intensity nodes [29] [see z-derivative plot in Fig. 4(a)]. To achieve a relatively constant localization precision along

*z*axis, a phase-shifting technique is used in such a way that the zero z-derivative of the PSF overlaps with non-zero z-derivative of the phase-shifted PSFs [16, 29]. However, this method requires multiple detectors [16]. Instead, our iMFM system simultaneously captures multiple focal shift tile images by using a single detector. The combined z-directive summed over all the tile images in iMFM [see Fig. 4(b)] allows us to achieve a relatively constant localization precision by using a single detector. When CA exists in the presence of 10nm emission spectrum, both MFM and iMFM localization precisions decrease. An iMFM polychromatic PSF can achieve average localization precision of (

*σ*) = (25.6nm, 24.0nm, 16.6nm) [cyan line], with 1.17-fold improvement laterally and 3.6-fold axially compared to MFM polychromatic PSF of localization precision (29.9nm, 28.3nm, 59.7nm) [green line]. Since the PSF is not symmetric along its

_{x}, σ_{y}, σ_{z}*x*and

*y*axis due to CA [see Fig. 3(d)], the localization precision are different in

*x*and

*y*directions. Those results indicate that iMFM could provide 3 to 4-fold improved axial localization precision in both aberrated and unaberrated systems.

#### 3.2. Initial point

In order to verify the theoretical analysis in section 3.1 and demonstrate the capability of iMFM to achieve higher axial localization precision than that of the single objective MFM, we have performed MLE reconstructions of single particle tracking using iMFM multifocal interferometric detection. However, it is known that the optimization problem in MLE as shown in Eq. (17) is non-convex, and a global minimum is not guaranteed to be found. Therefore, multiple random initial values are used for the MLE localization algorithm, and the optimal solution with the minimal cost function value is picked as the final reconstruction.

A better initial value closer to the global minimum than random initial values could improve the success and convergence of the MLE localization algorithm, but is difficult to be found since there is no prior information about the 3D position of the point in conventional microscopy imaging. In addition, it is impossible to tell whether the point is in the positive or negative defocus because the PSF is symmetric with respect to the focal plane and has the same blur size for the equal magnitude but opposite defocus.

Here, we propose a method to determine an initial axial position of a single point that is imaged by the MFM and iMFM microscopes. In MFM and iMFM, the PSF is not symmetric with respect to the focal plane since the points in the positive and negative defocus are focused in the opposite diffraction orders. Furthermore, each *z* position point is focused in different tiles, formulated by *z* = (*m* + *Bn*)Δ*z*, where (*m, n*) is the focused tile diffraction order, *B* and Δ*z* are the pre-designed parameters which are known priori. Therefore, if we can determine which tile image is most in focus by comparing the blur sizes of the tile-PSFs, then the initial axial position of the molecule *z*_{0} can be found as *z*_{0} = (*m* + *Bn*)Δ*z*. Also, since the focal step between two consecutive tiles is Δ*z*, the error distance between the initial estimation and the ground truth should be smaller than Δ*z*, i.e., $\left|{z}_{0}-\widehat{z}\right|\le \mathrm{\Delta}z$. In the simulation, we set $\left|{z}_{0}-\widehat{z}\right|\le 2\mathrm{\Delta}z$ in order to preserve both accuracy and speed of the MLE localization algorithm. To avoid ambiguities caused from points on the edge of a tile, a beam stop slightly smaller than the tile width can be placed at the intermediate image plane of both detection paths to restrict the lateral FOV [10]. An example of determining the initial axial position of a single point is shown in Fig. 6.

#### 3.3. MLE reconstruction

For each position $({\widehat{x}}_{0},{\widehat{y}}_{0},{\widehat{z}}_{0})$ of the molecule emitter, the acquired pixelated image under Poisson noise corruption is generated as

*C*denotes the pixel area on the detector plane, ${\widehat{h}}_{\text{iMFM}}$ is the normalized iMFM PSF with its integral equal to one, and

_{p}*N*

_{tot}is the total number of the photons collected by the iMFM PSF. The microscope and MFG parameters were the same as those used for PSF simulation in section 2.3 and CRLB calculation in section 3.1. The acquired image

*I*consisted of 3 × 3 tile images with a focal shift of Δ

*z*= 0.25

*µ*m. We assumed that each tile image had a region of interest (ROI) of 41 × 41 pixels, with pixel size of 4

*µ*m ×4

*µ*m. In addition, each pixel was composed of 4 × 4 sub-pixels for the purpose of the integral over the pixel area

*C*.

_{p}We simulated 50 images for each *z* position of the emitter between −1*µ*m and 1*µ*m from Eq. (21) and used MLE to back-calculate 3D position of the emitter, starting with a proper initial value as described in section 3.2. For each 3D position, a cluster of the positions containing 50 points was recovered and the mean squared error between the estimated position and true position was calculated. The estimation errors in the axial dimension by both MFM and iMFM are shown in Fig. 5(right) (black line and purple line, respectively). The simulation results indicate that the MLE estimation errors are well consistent with theoretical predictions (blue line and red line, respectively).

We also simulate a trajectory of a single emitter which follows a random walk from axial −1*µ*m toward 1*µ*m, shown in Fig. 7(a) in 3D view along with its projections onto the *xy, xz* and *yz* planes. The MFM and iMFM reconstructed trajectories are shown in Figs. 7(b) and 7(c), respectively. The difference between the ground truth and reconstructions along the lateral and axial directions is also plotted as a histogram of reconstruction residuals, shown in Figs. 7(d) and 7(e). The standard deviation for MFM axial resolution is 48.37nm and standard deviation for iMFM is 12.12nm, resulting in a 4-fold improvement in axial localization precision that iMFM detection PSF can offer over a large volume.

There are other two localization methods for multi-plane detection: (1) least-square fitting for bi-plane imaging and (2) 3D PSF Gaussian fitting for MFM. The MLE estimator used in the paper is more similar to the least-square fitting by assuming Poisson noise instead of Gaussian noise due to the limited-photon budge of the iMFM. The more details about the comparison between these methods can be found in [35, 36]. Another point worth mentioning here is that even though we only show single emitter tracking in the section, the proposed iMFM system is also applicable to multiple emitters tracking. In multiple emitters localization, we can first identify the emitters and then perform MLE routine for every identified particle, as introduced in least-squares fitting localization algorithm for bi-plane imaging [35].

#### 3.4. Optimal parameters design

So far our analysis and simulations assume a 10nm bandpass filter used in the chromatic aberrated iMFM system for the purpose of mitigation of CA. However, this narrow bandpass filter may not be optimal because it results in fewer photon detection and therefore decreases the localization precision for the iMFM. Here we plot the lateral and axial theoretical localization precision by assuming different bandwidth filters for both chromatic corrected iMFM [Fig. 8(a)] and chromatic aberrated iMFM system [Fig. 8(b)]. We consider four different bandwidths: 10nm (cyan), 20nm (red), 40nm (blue), and 80nm (black). For simplicity, the total number of the detected signal photon is assumed to be proportional to the bandwidth of the filter, i.e., 2500 signal photon for 10nm, 5000 for 20nm, 10000 for 40nm, and 20000 for 80nm. The background photon is 10 per pixel for all the cases. The other simulation parameters are same as those used in Fig. 5. The wavelength-dependent diffraction efficiency and both lateral and axial chromatic dispersion are also considered in the simulation. The average localization over the axial imaging range is also plotted as a function of the bandwidth for both iMFM systems in Fig. 8(c). From Fig. 8(c), we can see that the localization precision improves with the bandwidth of the filter even for the chromatically aberrated iMFM due to the increased total number of the collected photon. Note that the localization improvement in real experiment may not be as high as that in the simulation, because the total number of the detected photon is not linearly increased with the bandwidth of the filter due to the non-uniform distribution of the emission spectrum over different wavelengths.

In addition, we also investigate the effect of the focal spacing on the localization precision of the iMFM. In Fig. 9, we show the lateral (left) and axial localization precision of the chromatic corrected iMFM system with three different focal plane spacings, i.e., 100nm (red), 250nm (blue) and 400nm (black). In the simulation, we assume a 10nm bandwidth filter. The total number of the signal photon is 2500 and the background photon is 10 per pixel for all the cases. From Fig. 9, we can see that the effective axial imaging range of the iMFM can be enlarged by increasing the focal plane spacing at the expense of the severer oscillation of both lateral and axial localizations, which is consistent with the observation in MFM [26].

## 4. Isotropic 3D resolution for iMFM wide field imaging

#### 4.1. Nyquist sampling for iMFM wide-field imaging

For a band-limited system, the Nyquist sampling rate has to be satisfied in order to avoid aliasing. In our dual-objective iMFM microscope, the lateral cut-off frequency is 2*N A*/*λ* and the axial cut-off frequency is 2*n*_{0}/*λ*. Therefore, the Nyquist sampling distance in the sample space has to be equal or less than Δ* _{xy}* =

*λ*/(4

*N A*) in the lateral direction and Δ

*=*

_{z}*λ*/(4

*n*

_{0}) in the axial direction.

For an iMFM system with *N A* = 1.27, index of the refraction *n*_{0} = 1.338 and *λ* = 620nm, we found that Δ* _{xy}* = 122nm and Δ

*= 116nm. In the microscope design, the lateral sampling distance can be met by choosing a detector with proper pixel size such that ${d}_{\text{pixel}}/\widehat{M}\le {\mathrm{\Delta}}_{xy}$, and the axial sampling distance can be met by designing the MFG with a focal step |Δ*

_{z}*z*| ≤ Δ

*.*

_{z}#### 4.2. Numerical reconstruction

To confirm that iMFM provides axial super resolution and hence 3D isotropic resolution in wide field fluorescence imaging, we also performed reconstruction of a 3D synthetic extended object. We used the same parameter values for the microscope: NA = 1.27 with index of refraction 1.338, and magnification $\widehat{M}=100$. In order to record more than 9 focal images in a single shot, we designed a new MFG that produces an array of 5 × 5 focal shift images. The focal step was designed to be Δ*z* = 100nm in order to satisfy the Nyquist–Shannon sampling condition along *z*. The sensor size is assumed to be 1024 × 1024 pixels, with 12.2*µ*m × 12.2*µ*m pixel size. The size of each focal shift tile image is about 205 × 205 pixels. To avoid lateral convolution artifacts at the boundary, the lateral field of view (FOV) of a 3D synthetic extended object is confined to a central region of 129 × 129 pixels. The 3D synthetic object resembling the structure of the microtubules from [37], and its 3D spectrum are shown in Fig. 10(a). The 3D synthetic object is of size 129 × 129 × 49 voxels and each voxel size is 120nm × 120nm × 50nm.

The single shot 2D measurements for iMFM and MFM were simulated using Eq. (13). To simulate a microscope with chromatic aberrations (CA), we considered the emission bandwidth of 10nm by using a 10nm filter. The maximum number of photons detected by the brightest pixel was 500, and the corresponding Poisson noise was added in each measurement. For the reconstruction, we used Richardson-Lucy algorithm with total variation (TV) regularization [27], as discussed in section 2.5. The optimal regularization parameters were found by exhaustive search and the algorithm was run to converge, defined when the difference between two consecutive values of the cost function is smaller than a predefined threshold. Figures 10(b)–10(d) show the reconstruction with the use of MFM monochromatic PSF, iMFM monochromatic and polychromatic PSFs in the presence of CA, respectively. The fourth and fifth rows of Fig. 10 show the linecuts and spectra comparisons, respectively, between the ground truth, MFM, and iMFM reconstructions. We can see that two axial finer features separated by around 375nm axially are blurred to a single feature in the MFM reconstruction [Fig. 10(b)], but are well resolved by both unaberrated and aberrated iMFM reconstructions [Figs. 10(c) and 10(d)]. Figure 10 clearly demonstrates the ability of iMFM to recover sample high axial spatial frequencies beyond the detection cut-off of the single lens MFM, and therefore achieve super resolution in the axial direction even for a system with chromatic aberrations.

## 5. Conclusions

We proposed an interferometric multifocus microscopy (iMFM) system, a combination of multi-focus microscopy (MFM) with two opposing objective lenses in an interferometric microscope design, which provides a higher axial resolution than that of the single objective MFM and hence isotropic 3D resolution in a single shot. We examined the possibility and challenge of combining I^{2}M and MFM. The proposed iMFM involves employing two MFGs of opposite focal steps placed at the Fourier planes inside the interferometer. In this configuration, the emitted light from the same point source can be directed into the same tile on the detector along both paths of the interferometer, and therefore self-interfere after passing through two lenses. The mathematical formulations of the iMFM monochromatic and polychromatic PSFs were derived and the image formation models were provided. The iMFM PSFs were simulated to show that this new iMFM configuration is capable of recording multiple focal shift interferometry in a single exposure without focal scanning, significantly speeding up the acquisition process of conventional dual objective interferometric detection.

We demonstrate the proposed iMFM microscope with two applications: (i) single molecule tracking and (ii) wide field 3D extended object imaging. We calculated the Fisher information matrix (FIM) and the Cramér-Rao lower bound (CRLB) of iMFM for both monochromatic and polychromatic PSFs. The iMFM PSF contains almost uniformly high combined differential information along the optic axis due to the steepened axial features of interferometry and the simultaneously multifocal detection scheme, proving that the iMFM PSF is more effective for encoding a single point position than MFM PSF. For 2500 detected photons per objective, a background of 10 photons per pixel, MFGs of 3 × 3 tiles with a focal step of 0.25*µ*m and a total efficiency of 67%, and a single camera detection, which are typical conditions and values in practice [8, 32–34], the theoretical localization precision of (16.6nm, 16.6nm, 11.2nm) and (25.6nm, 24.0nm, 16.6nm) in three dimensions can be achieved for iMFM monochromatic and polychromatic PSFs, with a 4-fold and 3.6-fold axial resolution improvement compared with its MFM counterparts. To prevent the light loss due to the BS, a second detector can be introduced and placed at the other output of the BS, where resolution better than 10nm in three dimensions is obtainable. Another advantage of MFM and iMFM is the focal shift between tiles, which is known as a prior and can be used to estimate the axial initial values of 3D position to improve the accuracy and convergence of the MLE localization algorithms. The reconstruction errors by MLE with proposed initial value estimation are consistent with theoretical predictions.

For 3D wide field imaging, the dual-objective iMFM tile-OTF features about 4-fold enlarged support of transferred spatial frequencies in the axial direction compared with single lens MFM. Therefore, the iMFM microscope provides improved axial and more isotropic 3D resolution in wide field 3D extended object imaging in a single shot without focal scanning. We also show that even with chromatic aberrations in the presence of 10nm bandwidth emission, the iMFM is still capable of recovering high axial spatial frequencies beyond the detection cut-off. The axial resolution can be further improved if the dual objectives are also used for illumination in addition to collection as in I^{5}M and 4Pi C type microscopes.

Just as in I^{2}M/I^{5}M/4Pi dual-objective microscopes, the improved axial resolution in iMFM comes at the cost of system complexity, because two opposing lenses and two MFGs are needed for iMFM multifocal interferometric detection. However, we anticipate that the complexity of experimental implementation of our dual-objective iMFM is marginally greater to that of previously implemented multi-phase iPALM microscopes [16,31] where three beam splitters and three detectors as well as two opposing objective lenses are utilized. Compared with iPLAM, an unavoidable limitation of the iMFM multi-plane detection is a limited photon budget for each focal plane, because of (1) the photon loss from the MFG and beam-splitter, (2) photon splitting between different image planes and (3) chromatic dispersion in the chromatically aberrated system. However, the proposed iMFM can achieve a larger localization depths by simultaneously imaging multiple different focal planes. In addition, the iMFM system uses only one detector instead of three. Finally, our iMFM system can also be used for 3D extended object recovery without focal scanning, significantly speeding up the acquisition process of conventional dual-objective single plane interferometric detection in iPALM. We believe that the proposed iMFM microscope will be a useful tool in 3D dynamic event imaging where both high temporal and spatial resolution are required. The code is available upon request.

## Appendix

In this appendix, we show the derivation details of the multifocus grating (MFG) equation as shown in Eq. (6). Let’s start with a normal periodic grating *g* with spacing *d _{x}* and

*d*in the

_{y}*x*and

*y*directions, respectively. Because the grating is periodic and continuous, the Fourier transform (FT) of it yields discrete and aperiodic spectrum as

*u*and

*v*are spatial frequencies in

*x*and

*y*directions, respectively, and

*u*

_{0}= 1/

*d*and

_{x}*v*

_{0}= 1/

*d*are the intervals between consecutive samples in the discrete spectrum of the grating. In Fourier optics,

_{y}*u*=

*x*/(

*fλ*) and

*v*=

*y*/(

*fλ*) under the paraxial approximation, where

*x*and

*y*are spatial coordinates in detection plane. Note that the paraxial approximation holds here because the focal length

*f*of the relay system lens L4 in Fig. 1(b) is much larger than the sensor size.

For MFG, the geometrical distortions Δ* _{x}* and Δ

*are introduced in the grating pattern in the*

_{y}*x*and

*y*directions, respectively. Therefore, according to the FT shift theorem, the FT of the distorted grating

*g*

_{1}can be written as

*mµ*

_{0}and

*nv*

_{0}indicate different discrete spatial frequencies in

*x*and

*y*directions, respectively. As suggested in [8], the geometrical distortions are set to be Δ

*=*

_{x}*d*

_{x}n_{0}Δ

*z*cos

*θ*/

*λ*and Δ

_{c}*=*

_{y}*B*Δ

*to create a proper refocus. Therefore, Eq. (23) can be rewritten as*

_{x}*k*= 2

*πn*

_{0}/

*λ*,

*z*,

_{m}*= (*

_{n}*m*+

*Bn*) Δ

*z*, and

*x*

_{0}=

*f λ*/

_{c}*d*and

_{x}*y*

_{0}=

*f λ*/

_{c}*d*are spatial intervals between consecutive diffraction orders in the

_{y}*x*and

*y*directions, respectively on the detector plane for the emission central wavelength

*λ*. Equation (7) of the opposite focal shift MFG can also be derived in the similar way by setting the distortion to be −Δ

_{c}*and −Δ*

_{x}*in the two directions.*

_{y}## Funding

Biological Systems Science Division, Office of Biological and Environmental Research, Office of Science, U.S. Dept. of Energy, under Contract DE-AC02-06CH11357; National Science Foundation (NSF) CAREER grant IIS-1453192; Office of Naval Research (ONR) award N00014-15-1-2735.

## References

**1. **M. G. L. Gustafsson, D. A. Agard, and J. W. Sedat, “Sevenfold improvement of axial resolution in 3D wide-field microscopy using two objective-lenses,” Proc. SPIE **2412**, 147–156 (1995). [CrossRef]

**2. **P. Prabhat, S. Ram, E. S. Ward, and R. J. Ober, “Simultaneous imaging of different focal planes in fluorescence microscopy for the study of cellular dynamics in three dimensions,” IEEE Trans. Nanobioscience **3**(4), 237–242 (2004). [CrossRef]

**3. **S. Ram, P. Prabhat, J. Chao, E. S. Ward, and R. J. Ober, “High accuracy 3D quantum dot tracking with multifocal plane microscopy for the study of fast intracellular dynamics in live cells,” Biophys. J. **95**, 6025–6043 (2008). [CrossRef] [PubMed]

**4. **S. Ram, D. Kim, R. J. Ober, and E. S. Ward, “3D single molecule tracking with multifocal plane microscopy reveals rapid intracellular transferrin transport at epithelial cell barriers,” Biophys. J. **103**, 1594–1603 (2012). [CrossRef] [PubMed]

**5. **P. Prabhat, Z. Gan, J. Chao, S. Ram, C. Vaccaro, S. Gibbons, R. J. Ober, and E. S. Ward, “Elucidation of intracellular recycling pathways leading to exocytosis of the Fc receptor, FcRn, by using multifocal plane microscopy,” Proc. Natl. Acad. Sci. U.S.A. **104**, 5889–5894 (2007). [CrossRef] [PubMed]

**6. **P. M. Blanchard and A. H. Greenaway, “Simultaneous multiplane imaging with a distorted diffraction grating,” Appl. Opt. **38**(32), 6692–6699 (1999). [CrossRef]

**7. **P. A. Dalgarno, H. I. C. Dalgarno, A. Putoud, R. Lambert, L. Paterson, D. C. Logan, D. P. Towers, R. J. Warburton, and A. H. Greenaway, “Multiplane imaging and three dimensional nanoscale particle tracking in biological microscopy,” Opt. Express **18**, 877–884 (2010). [CrossRef] [PubMed]

**8. **S. Abrahamsson, J. Chen, B. Hajj, S. Stallinga, A. Y. Katsov, J. Wisniewski, G. Mizuguchi, P. Soule, F. Mueller, C. D. Darzacq, X. Darzacq, C. Wu, C. I. Bargmann, D. A. Agard, M. Dahan, and M. G. L. Gustafsson, “Fast multicolor 3D imaging using aberration-corrected multifocus microscopy,” Nat. Methods **10**, 60–63 (2013). [CrossRef]

**9. **S. Abrahamsson, M. McQuilken, S. B. Mehta, A. Verma, J. Larsch, R. Ilic, R. Heintzmann, C. I. Bargmann, A. S. Gladfelter, and R. Oldenbourg, “MultiFocus Polarization Microscope (MF-PolScope) for 3D polarization imaging of up to 25 focal planes simultaneously,” Opt. Express **23**, 7734–7754 (2015). [CrossRef] [PubMed]

**10. **X. Huang, A. Selewa, X. Wang, M. K. Daddysman, I. Gdor, R. Wilton, K. M. Kemner, S. Yoo, A. K. Katsaggelos, K. He, O. Cossairt, N. J. Ferrier, M. Hereld, and N. F. Scherer, “3D snapshot microscopy of extended objects,” https://arXiv:1802.01565 (2018).

**11. **J. Chen, Z. Zhang, L. Li, B. Chen, A. Revyakin, B. Hajj, W. Legant, M. Dahan, T. Lionnet, E. Betzig, R. Tjian, and Z. Liu, “Single-molecule dynamics of enhanceosome assembly in embryonic stem cells,” Cell , **156**(6), 1274–1285 (2014). [CrossRef] [PubMed]

**12. **J. Wisniewski, B. Hajj, J. Chen, G. Mizuguchi, H. Xiao, D. Wei, M. Dahan, and C. Wu, “Imaging the fate of histone Cse4 reveals de novo replacement in S phase and subsequent stable residence at centromeres,” eLife **3**, e02203 (2014). [CrossRef] [PubMed]

**13. **B. Hajj, J. Wisniewski, M. El Beheiry, J. Chen, A. Revyakin, C. Wu, and M. Dahan, “Whole-cell, multicolor superresolution imaging using volumetric multifocus microscopy,” Proc. Natl. Acad. Sci. U.S.A. **111**(49), 17480–17485 (2014). [CrossRef] [PubMed]

**14. **S. R. P. Pavani, M. A. Thompson, J. S. Biteen, S. J. Lord, N. Liu, R. J. Twieg, R. Piestun, and W. E. Moerner, “Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function,” Proc. Natl. Acad. Sci. U.S.A. **106**(9), 2995–2999 (2009). [CrossRef] [PubMed]

**15. **B. Bailey, D. L. Farkas, D. L. Taylor, and F. Lanni, “Enhancement of axial resolution in fluorescence microscopy by standing-wave excitation,” Nature **366**, 44–48 (1993). [CrossRef] [PubMed]

**16. **G. Shtengel, J. A. Galbraith, C. G. Galbraith, J. Lippincott-Schwartz, J. M. Gillette, S. Manley, R. Sougrat, C. M. Waterman, P. Kanchanawong, M. W. Davidson, R. D. Fetter, and H. F. Hess, “Interferometric fluorescent super-resolution microscopy resolves 3D cellular ultrastructure,” Proc. Natl. Acad. Sci. U.S.A. **106**(9), 3125–3130 (2009). [CrossRef] [PubMed]

**17. **D. Aquino, A. Schonle, C. Geisler, C. V. Middendorff, C. A. Wurm, Y. Okamura, T. Lang, S. W. Hell, and A. Egner, “Two-color nanoscopy of three-dimensional volumes by 4Pi detection of stochastically switched fluorophores,” Nat. Methods **8**(4), 353–359 (2011). [CrossRef] [PubMed]

**18. **I. Gdor, X. Wang, M. Daddysman, Y. Yifat, R. Wilton, M. Hereld, M. F. Noirot-Gros, and N. F. Scherer, “Particle Tracking by Repetitive Phase-Shift Interferometric Super Resolution Microscopy,” Opt. Lett. **43**, 2819–2822 (2018). [CrossRef] [PubMed]

**19. **M. Gu, *Advanced Optical Imaging Theory* (Springer, 1999, Chap. 6).

**20. **J. Bewersdorf, R. Schmidt, and S. W. Hell, “Comparison of I5M and 4Pi-microscopy,” J. Microsc. **222**(2), 105–117 (2006). [CrossRef] [PubMed]

**21. **M. Nagorni and S. W. Hell, “Coherent use of opposing lenses for axial resolution increase in fluorescence microscopy. I. Comparative study of concepts,” J. Opt. Soc. Am. A **18**(1), 36–48 (2001). [CrossRef]

**22. **W. Wang and G. Situ, “Interferometric rotating point spread function and its application in localization based super resolution imaging,” Sci. Rep. **7**(1), 5882 (2017). [CrossRef] [PubMed]

**23. **R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik **35**, 237–246 (1972).

**24. **J. N. Mait, “Understanding diffractive optic design in the scalar domain,” J. Opt. Soc. Am. A **12**, 2145–2158 (1995). [CrossRef]

**25. **B. Hajj, L. Oudjedi, J. B. Fiche, M. Dahan, and M. Nollmann, “Highly efficient multicolor multifocus microscopy by optimal design of diffraction binary gratings,” Sci. Rep. **12**(7), 5284 (2017). [CrossRef]

**26. **A. Tahmasbi, S. Ram, J. Chao, A. V. Abraham, F. W. Tang, E. S. Ward, and R. J. Ober, “Designing the focal plane spacing for multifocal plane microscopy,” Opt. Express **22**(14), 16706–16721 (2014). [CrossRef] [PubMed]

**27. **N. Dey, L. Blanc-Feraud, C. Zimmer, P. Roux, Z. Kam, J. C. Olivo-Marin, and J. Zerubia, “Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution,” Microsc. Res. Tech. **69**(4), 260–266 (2006). [CrossRef] [PubMed]

**28. **S. Yoo, P. Ruiz, X. Huang, K. He, M. Hereld, A. Selewa, M. Daddysman, N. Scherer, O. Cossairt, and A. K. Katsaggelos, “3D image reconstruction from multi-focus microscope: axial super-resolution and multiple-frame processing,” Proc. IEEE Int. Conf. Acoust. Speech, and Signal Processing, (2018).

**29. **C. V. Middendorff, A. Egner, C. Geisler, S. W. Hell, and A. Schönle, “Isotropic 3D Nanoscopy based on single emitter switching,” Opt. Express **16**(25), 20774–20788 (2008). [CrossRef]

**30. **S. Ram, P. Prabhat, E. S. Ward, and R. J. Ober, “Improved single particle localization accuracy with dual objective multifocal plane microscopy,” Opt. Express **17**(8), 6881–6898 (2009). [CrossRef] [PubMed]

**31. **S. W. Hell, R. Schmidt, and A. Egner, “Diffraction-unlimited three-dimensional optical nanoscopy with opposing lenses,” Nat. Photonics **3**, 381–387 (2009). [CrossRef]

**32. **R.J. Ober, S. Ram, and S.E. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. **86**, 1185–1200 (2004). [CrossRef] [PubMed]

**33. **T. Kues, R. Peters, and U. Kubitscheck, “Visualization and tracking of single protein molecules in the cell nucleus,” Biophys. J. **80**, 2954–2967 (2001). [CrossRef] [PubMed]

**34. **U. Kubitscheck, O. Kuckmann, T. Kues, and R. Peters, “Imaging and tracking single GFP molecules in solution,” Biophys. J. **78**, 2170–2179 (2000). [CrossRef] [PubMed]

**35. **M. J. Mlodzianoski, M. F. Juette, G. L. Beane, and J. Bewersdorf, “Experimental characterization of 3D localization techniques for particle-tracking and super-resolution microscopy,” Opt. Express **17**, 8264–8277 (2009). [CrossRef] [PubMed]

**36. **F. Aguet, D. V. D. Ville, and M. Unser, “A maximum-likelihood formalism for sub-resolution axial localization of fluorescent nanoparticles,” Opt. Express **13**, 10503–10522 (2005). [CrossRef] [PubMed]

**37. **D. Sage, L. Donati, F. Soulez, D. Fortun, G. Schmit, A. Seitz, R. Guiet, C. Vonesch, and M. Unser, “DeconvolutionLab2: An Open-Source Software for Deconvolution Microscopy,” Methods **115**, 28–41 (2017). [CrossRef] [PubMed]