Multiplexed phase-space imaging for 3D fluorescence microscopy

TL;DR

This paper introduces a multiplexed phase-space imaging technique for 3D fluorescence microscopy that enhances resolution, speed, and light efficiency over traditional scanning methods by using a spatial light modulator to pattern multiplexed apertures.

Contribution

The authors develop a novel multiplexed approach with a spatial light modulator to improve 3D fluorescence imaging, achieving diffraction-limited resolution and faster acquisition.

Findings

Reconstructed a 101-Megavoxel 3D volume with high resolution.

Demonstrated improved acquisition time and light throughput.

Enabled exploitation of sample sparsity for reduced data capture.

Abstract

Optical phase-space functions describe spatial and angular information simultaneously; examples include light fields in ray optics and Wigner functions in wave optics. Measurement of phase-space enables digital refocusing, aberration removal and 3D reconstruction. High-resolution capture of 4D phase-space datasets is, however, challenging. Previous scanning approaches are slow, light inefficient and do not achieve diffraction-limited resolution. Here, we propose a multiplexed method that solves these problems. We use a spatial light modulator (SLM) in the pupil plane of a microscope in order to sequentially pattern multiplexed coded apertures while capturing images in real space. Then, we reconstruct the 3D fluorescence distribution of our sample by solving an inverse problem via regularized least squares with a proximal accelerated gradient descent solver. We experimentally reconstruct a 101-Megavoxel 3D volume ($1010\times510\times500\mu$m with NA 0.4), demonstrating improved acquisition time, light throughput and resolution compared to scanning aperture methods. Our flexible patterning scheme further allows sparsity in the sample to be exploited for reduced data capture.

Equations

$$W(\mathbf{u},\mathbf{r})$$

$$\mathrm{MI}(\mathbf{u},\Delta \mathbf{u})\triangleq ⟨{\tilde{E}}^{\ast }(\mathbf{u}-\frac{\Delta \mathbf{u}}{2})\tilde{E}(\mathbf{u}+\frac{\Delta \mathbf{u}}{2})⟩,$$

$${\tilde{E}}_n({\mathbf{u}}_1;\alpha )=M_n(\lambda f{\mathbf{u}}_1){\tilde{E}}_s({\mathbf{u}}_1;\alpha ),$$

$$I_n(\mathbf{r};\alpha )$$

$$=\text{iiiint}{\mathrm{MI}}_n(\mathbf{u},\Delta \mathbf{u};\alpha )e^{i2\pi \Delta \mathbf{u}\cdot \mathbf{r}}{\mathrm{d}}^2\Delta \mathbf{u}{\mathrm{d}}^2\mathbf{u}$$

$$=\iint W_n(\mathbf{u},\mathbf{r};\alpha ){\mathrm{d}}^2\mathbf{u},$$

$${\tilde{I}}_n(\Delta \mathbf{u};\alpha )$$

$$=\iint M_n^{\ast }\left(\lambda f(\mathbf{u}-\frac{\Delta \mathbf{u}}{2})\right)M_n\left(\lambda f(\mathbf{u}+\frac{\Delta \mathbf{u}}{2})\right){\mathrm{MI}}_s(\mathbf{u},\Delta {\mathbf{u}}^{\prime };\alpha ){\mathrm{d}}^2\mathbf{u},$$

$$I_n(\mathbf{r})$$

$$\tilde{E}_{s}({\mathbf{u}}_{1};{\mathbf{r}}_{s},z_{s},\lambda)=\begin{cases}e^{i\frac{2\pi}{\lambda}(-z_{s})\sqrt{1-\lambda^{2}|{\mathbf{u}}_{1}|^{2}}-i2\pi{\mathbf{r}}_{s}\cdot{\mathbf{u}}_{1}},&\mbox{$|{\mathbf{u}}_{1}|<\frac{NA}{\lambda}$}\\ 0,&\mbox{otherwise}\end{cases}.$$

$$I_n(\mathbf{r};{\mathbf{r}}_s,z_s,\lambda )=\iint K_{M_n,z_s}(\Delta \mathbf{u};\lambda )e^{i2\pi (\mathbf{r}-{\mathbf{r}}_s)\cdot \Delta \mathbf{u}}{\mathrm{d}}^2\Delta \mathbf{u},$$

$$\begin{array}{rl}{K}_{M_n,z_s}(\Delta \mathbf{u};\lambda )=\iint & M_n^{\ast }\left(\lambda f(\mathbf{u}-\frac{\Delta \mathbf{u}}{2})\right)M_n\left(\lambda f(\mathbf{u}+\frac{\Delta \mathbf{u}}{2})\right)\\ & e^{i\frac{2\pi z_s}{\lambda }\left(\sqrt{1-{\lambda }^2|\mathbf{u}-\frac{\Delta \mathbf{u}}{2}{|}^2}-\sqrt{1-{\lambda }^2|\mathbf{u}+\frac{\Delta \mathbf{u}}{2}{|}^2}\right)}{\mathrm{d}}^2\mathbf{u}\end{array}$$

$$I_n(\mathbf{r})=\iiint \left(\iint \sum _{\lambda }S(\lambda )K_{M_n,z_s}(\Delta \mathbf{u};\lambda )e^{i2\pi (\mathbf{r}-{\mathbf{r}}_s)\cdot \Delta \mathbf{u}}{\mathrm{d}}^2\Delta \mathbf{u}\right)C({\mathbf{r}}_s,z_s){\mathrm{d}}^2{\mathbf{r}}_s\mathrm{d}{z}_s.$$

$$\mathbf{y}=\mathbf{A}\mathbf{c}.$$

$${\mathrm{argmin}}_{\mathbf{c}\ge 0}\frac{1}{2}\Vert \mathbf{y}-{\mathbf{y}}_{\mathrm{meas}}{\Vert }_2^2+\mu \Vert \mathbf{W}\mathbf{c}{\Vert }_1,$$

Full text

OCIS codes: (100.3010) Image reconstruction techniques, (110.0180) Microscopy, (110.4980) Partial coherence in imaging, (110.1758) Computational imaging

References and links

• [1]

I. Cox, C. Sheppard, and T. Wilson, “Super-resolution by confocal fluorescent microscopy,” Optik 60, 391–396 (1982).

• [2]

F. Helmchen and W. Denk, “Deep tissue two-photon microscopy,” Nat. Methods 2, 932–940 (2005).

• [3]

T. Planchon, L. Gao, D. Milkie, M. Davidson, J. Galbraith, C. Galbraith, and E. Betzig, “Rapid three-dimensional isotropic imaging of living cells using bessel beam plane illumination,” Nat. Methods 8, 417–423 (2011).

• [4]

M. Levoy, R. Ng, A. Adams, M. Footer, and M. Horowitz, “Light field microscopy,” in “SIGGRAPH ’06,” (New York, NY, USA, 2006), pp. 924–934.

• [5]

M. Broxton, L. Grosenick, S. Yang, N. Cohen, A. Andalman, K. Deisseroth, and M. Levoy, “Wave optics theory and 3-d deconvolution for the light field microscope,” Opt. Express 21, 25418–25439 (2013).

• [6]

Y. Y. Schechner, R. Piestun, and J. Shamir, “Wave propagation with rotating intensity distributions,” Phys. Rev. E 54, R50–R53 (1996).

• [7]

S. R. P. Pavani, M. A. Thompson, J. S. Biteen, S. J. Lord, N. Liu, R. J. Twieg, R. Piestun, and W. Moerner, “Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix point spread function,” Proceedings of the National Academy of Sciences 106, 2995–2999 (2009).

• [8]

Y. Shechtman, S. J. Sahl, A. S. Backer, and W. E. Moerner, “Optimal point spread function design for 3D imaging,” Phys. Rev. Lett. 113, 133902 (2014).

• [9]

N. Ji, D. E. Milkie, and E. Betzig, “Adaptive optics via pupil segmentation for high-resolution imaging in biological tissues,” Nat. Methods 7, 141–147 (2010).

• [10]

M. Testorf, B. Hennelly, and J. Ojeda-Castañeda, Phase-Space Optics: Fundamentals and Applications (McGraw-Hill Professional, 2009).

• [11]

M. Bastiaans, “The Wigner distribution function applied to optical signals and systems,” Opt. Commun. 25, 26 – 30 (1978).

• [12]

A. Accardi and G. Wornell, “Quasi light fields: extending the light field to coherent radiation,” J. Opt. Soc. Am. A 26, 2055–2066 (2009).

• [13]

Z. Zhang and M. Levoy, “Wigner distributions and how they relate to the light field,” in “Computational Photography (ICCP), 2009 IEEE International Conference on,” (2009), pp. 1–10.

• [14]

R. Ng, M. Levoy, M. Bredif, G. Duval, M. Horowitz, and P. Hanrahan, “Light field photography with a hand-held plenoptic camera,” Tech. Rep. CTSR 2005-02, Stanford (2005).

• [15]

R. Ng and P. Hanrahan, “Digital correction of lens aberrations in light field photography,” in “Contract Proceedings 2006,” (International Society for Optics and Photonics, 2007), pp. 63421E–63421E.

• [16]

N. C. Pégard, H.-Y. Liu, N. Antipa, M. Gerlock, H. Adesnik, and L. Waller, “Compressive light-field microscopy for 3D neural activity recording,” Optica 3, 517–524 (2016).

• [17]

D. L. Marks, R. A. Stack, and D. J. Brady, “Three-dimensional coherence imaging in the Fresnel domain,” Appl. Opt. 38, 1332–1342 (1999).

• [18]

J. K. Wood, K. A. Sharma, S. Cho, T. G. Brown, and M. A. Alonso, “Using shadows to measure spatial coherence,” Opt. Lett. 39, 4927–4930 (2014).

• [19]

D. Christodoulides, E. Eugenieva, T. Coskun, M. Segev, and M. Mitchell, “Equivalence of three approaches describing partially incoherent wave propagation in inertial nonlinear media,” Phys. Rev. E 63, 35601 (2001).

• [20]

L. Tian, Z. Zhang, J. C. Petruccelli, and G. Barbastathis, “Wigner function measurement using a lenslet array,” Opt. Express 21, 10511–10525 (2013).

• [21]

L. Waller, G. Situ, and J. Fleischer, “Phase-space measurement and coherence synthesis of optical beams,” Nat. Photonics 6, 474–479 (2012).

• [22]

H.-Y. Liu, E. Jonas, L. Tian, J. Zhong, B. Recht, and L. Waller, “3D imaging in volumetric scattering media using phase-space measurements,” Opt. Express 23, 14461–14471 (2015).

• [23]

S. Cho, M. A. Alonso, and T. G. Brown, “Measurement of spatial coherence through diffraction from a transparent mask with a phase discontinuity,” Opt. Lett. 37, 2724–2726 (2012).

• [24]

H.-Y. Liu, J. Zhong, and L. Waller, “4D phase-space multiplexing for fluorescent microscopy,” in “SPIE BiOS,” (International Society for Optics and Photonics, 2016), pp. 97200A–97200A.

• [25]

A. Levin, R. Fergus, F. Durand, and W. T. Freeman, “Image and depth from a conventional camera with a coded aperture,” ACM transactions on graphics (TOG) 26, 70 (2007).

• [26]

C. Liang, T. Lin, B. Wong, C. Liu, and H. Chen, “Programmable aperture photography: Multiplexed light field acquisition,” in “ACM SIGGRAPH 2008 Papers,” (2008), pp. 55:1–55:10.

• [27]

P. Green, W. Sun, W. Matusik, and F. Durand, “Multi-aperture photography,” in “ACM Transactions on Graphics (TOG),” , vol. 26 (ACM, 2007), vol. 26, p. 68.

• [28]

J. Chang, I. Kauvar, X. Hu, and G. Wetzstein, “Variable aperture light field photography: Overcoming the diffraction-limited spatio-angular resolution tradeoff,” in “Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition,” (2016), pp. 3737–3745.

• [29]

A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM Journal on Imaging Sciences 2, 183–202 (2009).

• [30]

A. Walther, “Radiometry and coherence,” J. Opt. Soc. Am. 58, 1256–1259 (1968).

• [31]

J. Goodman, Introduction to Fourier optics (Roberts & Co., 2005).

• [32]

M. Bertero and P. Boccacci, “A simple method for the reduction of boundary effects in the Richardson-Lucy approach to image deconvolution,” Astronomy and Astrophysics 437, 369–374 (2005).

1 Introduction

3D fluorescence microscopy is a critical tool for bioimaging, since most samples are thick and can be functionally labeled. High-resolution 3D imaging typically uses confocal [1], two-photon [2] or light sheet microscopy [3]. Because these methods all involve scanning, they are are inherently limited in terms of speed or volume. Light field microscopy [4, 5], on the other hand, achieves single-shot 3D capture, but sacrifices resolution. High resolution and single-shot capture are possible with coded aperture microscopy [6, 7, 8, 9]; however, this requires an extremely sparse sample. Here, we describe a multi-shot coded aperture microscopy method for high-resolution imaging of large and dense volumes with efficient data capture.

Our work fits into the framework of phase-space optics, which is the general term for any space-angle description of light [10, 11, 12, 13]. Light fields, for example, are phase-space functions for ray optics with incoherent light. The light field describes each ray’s position and angle at a particular plane, which can be used for digital refocusing [14], 3D imaging, aberration correction [15] or imaging in scattering [16]. In microscopy, wave-optical effects become prominent and so ray optics no longer suffices if one wishes to achieve diffraction-limited resolution. Hence, we use Wigner functions, which are the wave-optical analog to light fields [13, 12].

The Wigner function describes spatial and spatial frequency (propagation angle) information for a wave-field of arbitrary coherence. It converges to the light field in the limit of incoherent ray optics [13]. Capturing Wigner functions is akin to spatial coherence imaging [17, 18], since they contain the same information as mutual intensity and coherent mode decompositions [10, 19]. In the case of fluorescence microscopy, where the object is a set of incoherent emitters, Wigner representations define the incoherent 3D Optical Transfer Function (OTF) [10].

The 4D nature of phase space (two spatial and two spatial-frequency dimensions) poses significant measurement challenges. Single-shot schemes (e.g. lenslet arrays [20]) must distribute the 4D function across a 2D sensor, severely restricting resolution. Scanning aperture methods [21, 22, 23] are slow ($\sim$minutes) and light inefficient. We seek here a flexible trade-off between capture time and resolution, with the ability to exploit sparsity for further data reduction.

Our experimental setup consists of a widefield fluorescence microscope with a spatial light modulator (SLM) in Fourier space (the pupil plane). The SLM implements a series of quasi-random coded aperture patterns, while collecting real space images for each (Fig. 1) [24]. The recovered 4D phase space has very large pixel count (the product of the pixel counts of the SLM and the sensor) $\sim 10^{12}$. Compared to scanning aperture methods [21, 22], the new scheme has three major benefits. First, it achieves better resolution by capturing high-frequency interference effects (high-order correlations). This enables diffraction-limited resolution at the microscope’s full numerical aperture (NA). Second, we achieve higher light throughput by opening up more of the pupil in each capture; this can be traded for shorter exposure time and faster acquisition. Third, the multiplexed nature of the measurements means that we can employ compressed sensing approaches (when samples are sparse) in order to capture fewer images without sacrificing resolution. This means that the number of images required scales not with the reconstructed number of resolved voxels, but rather with the sparsity of the volume.

Our method can be thought of as a multi-shot coded aperture scheme for diffraction-limited 3D fluorescence microscopy. It is analogous to coded aperture photography [25, 26, 27, 28]; however, we use a wave-optical model to account for diffraction effects, so intensity measurements are nonlinear with complex-field. Fluorescent imaging allows a simplification of the forward model, since each fluorophore is spatially coherent with itself but incoherent with all other fluorophores. Our reconstruction algorithm then becomes a large-scale inverse problem akin to multi-image 3D deconvolution, formulated as a convex $\ell_{1}$ regularized least-square error problem and solved by a fast iterative shrinkage-thresholding algorithm (FISTA) [29].

2 Theory and methods

We use the Wigner function (WF) and its close relative, the Mutual Intensity (MI) function, to describe our multiplexing scheme. The WF describes position and spatial frequency, where spatial frequency can be thought of as the direction of ray propagation in geometrical optics. The concept was introduced by Walther [30] to describe the connection between ray intensity and wave theory. Here we assume that the light is quasi-monochromatic, temporally stationary and ergodic, so the Wigner function is [11, 10]:

[TABLE]

where ${\mathbf{r}}=(x,y)$ denotes transverse spatial coordinates, ${\mathbf{u}}=(u_{x},u_{y})$ denotes spatial frequency coordinates and $\langle\cdot\rangle$ denotes ensemble average. The quantity contained by angle brackets,

[TABLE]

is, apart from a coordinate transform, the Mutual Intensity. $E(\mathbf{r})$ is a spatially coherent electric field (e.g. from a single fluorophore) and $\tilde{E}({\mathbf{u}})$ is its Fourier transform. The ensemble average allows representation of both coherent and partially (spatially) coherent light. Here, we assume that the object is a 3D volume of incoherent emitters with no occlusions. Thus, the phase-space description of light from the object is a linear sum of that from each fluorophore.

2.1 General forward model

Our forward model relates each coded aperture’s captured image to the 3D object’s intensity. Each image contains information from multiple frequencies and their interference terms (which are key to resolution enhancement). The MI framework facilitates our mathematical analysis, since the projection of a mask in MI space is equivalent to applying the 3D Optical Transfer Function (OTF) to an incoherent source (see Fig. 2). MI analysis further reveals the interference effects that may be obscured by looking only at a projected quantity (the OTF). For simplicity, we first describe the forward model of a point source and then generalize to multiple mutually incoherent sources.

Consider a complex electric field, specified by some properties, $\alpha$ (e.g. location and wavelength of a point source). The field $E_{s}({\mathbf{r}}_{1};\alpha)$ acts like a unique coherent mode – it interferes coherently with itself but not with other modes. For a single coherent mode at the front focal plane of the $4f$ system in Fig. 1, the complex-field just before the SLM is the Fourier transform of the input complex-field [31], expressed as $\tilde{E}_{s}({\mathbf{u}}_{1},\alpha)$, where ${\mathbf{u}}_{1}$ is in units of spatial frequency that relate to the lateral coordinates ${\mathbf{x}}_{1}$ on the SLM by multiplying the wavelength $\lambda$ and the front focal length $f$, such that ${\mathbf{x}}_{1}=\lambda f{\mathbf{u}}_{1}$. The field $\tilde{E}_{s}({\mathbf{u}}_{1};\alpha)$ is then multiplied by the coded aperture, giving,

[TABLE]

where $\tilde{E}_{n}$ is the patterned complex-field in Fourier space and $M_{n}$ represents the $n$’th binary coded aperture pattern. At the sensor plane (the back focal plane of the $4f$ system) the intensity $I_{n}$ is:

[TABLE]

where ${\mathbf{u}}_{2}$ is a duplicate of ${\mathbf{u}}_{1}$. The intensity image can be alternately related to the MI or WF by a simple coordinate transform ${\mathbf{u}}_{1}={\mathbf{u}}+\Delta{\mathbf{u}}/2$, ${\mathbf{u}}_{2}={\mathbf{u}}-\Delta{\mathbf{u}}/2$:

[TABLE]

where $\mathrm{MI}_{n}$ and $W_{n}$ are the MI and WF associated with the field modified by the $n$’th mask.

Describing the intensity measurement in terms of both its WF and MI gives new insights. Eq. (7) shows that the intensity image is a projection of the patterned Wigner function $W_{n}$ across all spatial frequencies. Alternately, looking at the Fourier transform of $I_{n}$, denoted by $\tilde{I}_{n}$, we can interpret the intensity measurement as a patterning and projection of the source MI:

[TABLE]

where $\mathrm{MI}_{s}$ is the MI of source field $E_{s}$. This interpretation is illustrated for a 1D complex-field in Fig. 2. Each SLM pattern probes different parts of the input MI. By using multiple coded apertures, one may reconstruct the MI from its projections.

The extension of our forward model to multiple coherent modes is straightforward. Since each mode is mutually incoherent with others, we sum over all of them with their weights $C(\alpha)$:

[TABLE]

2.2 Forward model for incoherent sources

From the general forward model of Eq. (9), we can now derive the case in which the object is a 3D incoherent intensity distribution (e.g. a fluorescent sample). We specify $\tilde{E}_{s}$ in Eq. (3) to be from a single-color fluorophore located at position $({\mathbf{r}}_{s},z_{s})$ where $z_{s}$ is the defocus distance. The parameter $\alpha$ is $({\mathbf{r}}_{s},z_{s},\lambda)$ for a single point source, and our goal is to solve for the coefficients $C(\alpha)$ which represents intensity of each. To account for off-focus point sources, the field can be propagated to $-z_{s}$ using angular spectrum propagation [31]:

[TABLE]

Using Eqs. (2), (6) and (8), after some algebra the intensity becomes:

[TABLE]

where

[TABLE]

is the kernel for mask $M_{n}$ at depth $z_{s}$. Plugging Eq. (11) into Eq. (9) gives the final expression of the forward model. Before doing so, we assume that the fluorescent color spectrum $S(\lambda)$ is identical and known for all flourophores. Thus, we can decompose the mode weights in Eq. (9) into a product of the spectral weight and the 3D intensity distribution, $S(\lambda)C({\mathbf{r}}_{s},z_{s})$, giving:

[TABLE]

Equation (13) describes the forward model for a 3D fluorescent object $C({\mathbf{r}}_{s},z_{s})$ with no occlusions. The term in parentheses is a convolution kernel describing the 3D Point Spread Function (PSF) for mask $M_{n}$ (shown in Fig. 1). For simplicity, we assume here no scattering, though incorporating the scattering forward model in [22, 16] is straightforward.

2.3 Inverse problem

Based on the raw data and forward model, the inverse problem is formulated as a nonlinear optimization. Our goal is to reconstruct the 3D intensity distribution $C({\mathbf{r}}_{s},z_{s})$ from the measured images. To do so, we aim to minimize data mismatch, with an $\ell_{1}$ regularizer to mitigate the effects of noise (and promote sparsity where applicable). The mismatch is defined as the least-square error between the measured intensity images and the intensity predicted by our forward model (Eq. (13)). This formulation has a smooth part and a non-smooth part in the objective function and is efficiently solved by a proximal gradient descent solver (FISTA [29]).

To formulate the inverse problem, we first discretize the forward model in Eq. (13) to be

[TABLE]

Here ${\bm{y}}\in\mathbb{R}^{MP\times 1}$ corresponds to predicted images on the sensor; each small chunk ($\in\mathbb{R}^{P\times 1}$) of ${\bm{y}}$ is a vectorized image $I_{n}({\mathbf{r}})$. We discretize ${\mathbf{r}}$ into $P$ pixels, and the number of masks is $M$, so $n=1\ldots M$. Similarly, we discretize ${\mathbf{r}}_{s}$ and $z_{s}$ into $P^{\prime}$ pixels and $L$ samples, respectively, to obtain a vectorized version of $C({\mathbf{r}}_{s},z_{s})$ (${\bm{c}}\in\mathbb{R}^{LP^{\prime}\times 1}$). The matrix ${\bm{A}}\in\mathbb{R}^{MP\times LP^{\prime}}$, which is not materialized, represents the summation and convolution in Eq. (13) using 2D Fast Fourier Transforms (FFTs) for each subvector ($\in\mathbb{R}^{P^{\prime}\times 1}$) of ${\bm{c}}$, with zero-padding to avoid periodic boundary condition errors. The convolution kernel is precomputed and stored for speed.

The inverse problem becomes a data fidelity term plus an $\ell_{1}$ regularization parameter $\mu$:

[TABLE]

where ${\bm{y}}_{\mathrm{meas}}\in\mathbb{R}^{MP\times 1}$ is the measured intensity. We also use a diagonal matrix ${\bm{W}}\in\mathbb{R}^{LP^{\prime}\times LP^{\prime}}$ to lower the weight of point sources near the borders of images whose light falls off the sensor. Each diagonal entry of ${\bm{W}}$ is obtained by summing the corresponding column in ${\bm{A}}$. Outside point sources may also contribute to the measured intensity due to defocus; hence, we use an extended field-of-view method [32] to solve for more sample points in ${\bm{c}}$ than ${\bm{y}}$ (i.e. $P^{\prime}>P$).

3 Design of coded apertures

In the scanning-aperture scheme [21, 22], smaller apertures give better frequency sampling of the 4D phase space, at a cost of: 1) lower resolution, 2) lower signal-to-noise ratio (SNR) and 3) large data sets. Our multiplexing scheme alleviates all of these problems. Multiplexing achieves diffraction-limited resolution by additionally capturing interference terms, which cover the full NA-limited bandwidth. This is evident in the Fourier transform of the captured images (Fig. 3). The SNR improvement is also visible; the multiplexed image is less noisy.

Our masks are chosen by quasi-random non-replacement selection. We section the SLM plane into 18$\times$18 square blocks and keep only the 240 blocks that are inside the system NA. For each mask, we open 12 blocks, selected randomly from the blocks remaining after excluding ones that were open in previous sequences. In this scheme, the full NA can be covered by 20 masks. To allow for both diversity and redundancy, we choose to to cover the entire pupil 5 times, resulting in 100 multiplexed aperture patterns, one of which is shown in Fig. 3.

Importantly, the number of multiplexed patterns can be flexibly chosen to trade off accuracy for speed of capture. For instance, by increasing the number of openings in each pattern, we can cover the entire pupil with fewer patterns. This means that we may be able to reconstruct the object from fewer measurements, if the inverse problem is solvable. By using a priori information about the object (such as sparsity in 3D) as a constraint, we can solve under-determined problems with fewer total measured pixels than voxels in the reconstruction.

4 Experiments

Our experimental setup consists of the $4f$ system ($f_{1}=250~{}mm$,$f_{2}=225~{}mm$) shown in Fig. 1, with an additional 4$f$ system in front, made of an objective lens (20$\times$ NA 0.4) and a tube lens ($f=200~{}mm$) to image the sample at the input plane. The SLM (1400$\times$1050 pixels of size $10.3{\mu\mathrm{m}}$) is a liquid crystal chip from a 3-LCOS projector (Canon SX50) which is reflective and polarization-sensitive, so we fold the optical train with a polarization beam splitter and insert linear polarizers. Our sensor (Hamamatsu ORCA-Flash4.0 V2) captures the multiplexed images and synchronizes with the SLM via computer control.

Our sample is a fixed fluorescent brine shrimp (Carolina Biological). It is relatively dense, yet does not fill the entire 3D volume. The reconstructed 3D intensity (Fig. 4(a), 4(b) and 4(f)-4(h)) is stitched from five volume reconstructions, each with 640$\times$640$\times$120 voxels to represent the sample volume of 455$\times$455$\times$600$~{}{\mu\mathrm{m}}$. The reconstruction is cropped to the central part of our extended field of view, so the final volume contains 1422$\times$715$\times$100 voxels corresponding to 1010$\times$510$\times$500${\mu\mathrm{m}}$. The dataset size is large (9 GB), and since the size of the 3D array is $\sim$5$\times 10^{7}$ without the extended field-of-view and the measured data is $\sim$4$\times 10^{7}$, the number of operations for evaluating Eq. (13) is on the order of 3$\times 10^{10}$. This takes 4 seconds to compute on a computer with 48-core 3.0 GHz CPUs and requires 94 GB memory to store the kernel (Eq. (12)).

The reconstructed 3D intensity is shown in Fig. 4, alongside images from a confocal microscope and a widefield focus stack, for comparison. Both our method and the focus stack use a 0.4 NA objective, while the confocal uses 0.25 NA; hence, the confocal results should have slightly better resolution. Our reconstructed slices appear to have slightly lower resolution than the defocus stack and confocal, possibly due to the missing information in the frequency mutual intensity illustrated in Fig. 2(c), which is greatly undersampled. As expected, the depth slices of our reconstruction have better rejection of information from other depths, similar to the confocal images.

To illustrate the flexible tradeoff between capture time (number of coded apertures used) and quality, we show reconstructions in Fig. 5 using different numbers of coded aperture images. The case of only 1 image corresponds to a single coded aperture and gives a poor result, since the sample is relatively dense. However, with as few as 10 images we obtain a reasonable result, despite the fact that we are solving a severely under-determined problem. This is possible because the measurements are multiplexed and so the $\ell_{1}$ regularizer acts as a sparsity promoter.

5 Conclusion

We demonstrated 3D reconstruction of a large-volume high-resolution dense fluorescent object from multiplexed phase-space measurements. An SLM in Fourier space dynamically implements quasi-random coded apertures while intensity images are collected in real space for each coded aperture. Theory is developed in the framework phase space, with relation to Mutual Intensity functions and 3D OTFs. Reconstruction is formulated as an $\ell_{1}$-regularized least-square problem. This method enables diffraction-limited 3D imaging with high resolution across large volumes, efficient data capture and a flexible acquisition scheme for different types and sizes of samples.

Funding

The Office of Naval Research (ONR) (Grant N00014-14-1-0083).

ACKNOWLEDGMENTS

The authors thank Eric Jonas, Ben Recht, Jingzhao Zhang and the AMP Lab at UC Berkeley for help with computational resources.