1. Introduction

3D quantitative phase imaging (QPI) has found use in many biological applications by providing 3D refractive index (RI) information of the samples [1]. Optical diffraction tomography (ODT) [2] is the most widely used technique for 3D QPI. It works by directly measuring the complex field using an interferometry setup from multiple angles and then recovering the 3D RI distribution with an inverse scattering model. However, the need for coherent illumination and an interferometric path complicates the system and makes it prone to speckle artifacts [3].

Intensity diffraction tomography (IDT) is an alternative 3D QPI technique that recovers 3D phase information from intensity-only measurements on a “non-interferometry” setup [4]. Recently, several groups have demonstrated that IDT can be easily implemented on a standard microscope using a programmable LED array [5–13]. Most of the IDT acquisition strategies require taking hundreds of images to achieve high-resolution 3D RI reconstruction. Recently, our group has developed two strategies to push the acquisition speed and enabled visualizing dynamic biological samples. The annular IDT achieved a 10 Hz volume rate by using an LED ring illuminator that provides illumination angles matching the objective’s numerical aperture (NA) [7]. The multiplexed IDT provided a $4 - {10}\;\textrm{Hz}$ acquisition rate by simultaneously tuning on multiple LEDs using the most widely used LED matrix [10].

The 3D reconstruction of IDT requires an inverse scattering model, similar to ODT. The single-scattering models are based on the first-Born [5,7,10] or first-Rytov [12] approximation, allowing closed-form solutions due to the linear relationship between the scattering potential and the measurements. However, the accuracy of these models are fundamentally limited by the weak scattering assumption. Recently, several multiple-scattering models have been developed based on the beam propagation method (BPM) [6,14] and multi-layer Born (MLB) model [8] and showed enhanced resolution and fidelity on strongly scattering biological samples, such as C. elegans embryos and worm. However, the accuracy of the BPM degrades for high-NA systems since it relies on the paraxial approximation [15], which limits its effectiveness for high-resolution imaging using high-NA optics. The MLB model relies on the local weak scattering approximation within each layer, which limits the model fidelity for strongly scattering samples. Recently, a split-step non-paraxial (SSNP) multiple-scattering model has been developed for ODT [16], which demonstrated substantial improvement in reconstruction quality compared to the single-scattering and BPM methods. This ODT setup, however, needs a laser and a spatial light modulator to provide the coherent illumination and a interferometric path to retrieve the phase information of the light field. Besides, the large number of measurements used limits this system’s ability of high-speed imaging for dynamic samples. Inspired by this work, here we develop an SSNP IDT model for the two fast-acquisition IDT modalities. We demonstrate that our model can effectively provide high-quality 3D RI reconstruction on several biological samples with simple “non-interferometric” IDT setups.

The SSNP model is a wave propagation model derived from the Helmholtz equation [17]. Similar to the BPM and MLB, the SSNP works by axially splitting the sample into multiple slices and propagate the wave slice-wise. The BPM and MLB rely on local homogeneous RI and weak RI contrast, respectively, to decouple the diffraction from the local phase modulation / single scattering within each slice. Instead, the SSNP simultaneously propagates both the field and the axial derivative of the field so that the coupling between the diffraction and local phase modulation is modeled jointly without making assumptions about the local RI distribution [16]. Due to the same “multi-slice” structure as the BPM and MLB, together with the fast Fourier transform (FFT)-based propagation implementation [16], the SSNP incurs similar computational costs while providing improved fidelity on high-scattering biological samples.

Here we first extend the SSNP model to IDT and apply it to two IDT modalities, including the sequential and multiplexed IDT (as summarized in Fig. 1(a)-(d). Next, we derive an efficient SSNP-IDT forward model (as illustrated in Fig. 1(e)). We then derive a reconstruction algorithm by solving a regularized least-squares optimization (Fig. 1(f)). We elucidate on how the analytical gradient of the SSNP-IDT model can be efficiently computed, akin to “back-propagation”. Furthermore, we devise a unified algorithmic framework that allows fast 3D recovery from both sequential and multiplexed IDT measurements based on a memory efficient automatic differentiation framework. To quantitatively evaluate the fidelity of the forward model and reconstruction algorithm, we conduct numerical simulations. Lastly, we experimentally demonstrate high-quality, large field-of-view (FOV) 3D RI reconstructions on buccal epithelial cells and a live C. elegans worm from annular IDT and live C. elegans embryos from multiplexed IDT.

 

Fig. 1. (a-b) Our annular IDT setup uses an LED ring of radius $r$ to illuminate the sample sequentially. The ring is placed $h$ away from the sample and the illumination angle $\theta =\arctan (r/h)$ matches the objective NA. (c-d) Our multiplexed IDT setup uses an LED matrix and illuminates the sample using several LEDs simultaneously. (e) The SSNP model involves successive propagation $\mathbf {P}$ and scattering $\mathbf {Q}$ operations to compute the scattered field and its axial derivative $\boldsymbol {\Phi }_i$ from one slice to the next. (f) 3D reconstruction by iteratively solving the SSNP-based optimization problem.

Download Full Size | PDF

2. Theory

2.1 SSNP theory

The SSNP model first discretizes the 3D sample into a series of axial ($z$) slices and then calculates the internal field slice-by-slice ($xy$), as illustrated in Fig. 1(e). In the following, we provide the key steps to derive the governing equations for light propagating through the sample by following [16,17], and then describe the integration to specific IDT systems in separate sections.

The scalar wave propagation in an inhomogeneous medium with a 3D RI distribution $n(\boldsymbol {r})$ satisfies the Helmholtz equation

Next, we describe the numerical method to efficiently calculate Eq. (2)–(4). First, we discretize the sample volume as a series of $z$-slices. Second, to compute the field propagating a small distance $\Delta z$ between two adjacent slices, we approximate Eq. (2) as a first-order homogeneous linear system of differential equations with constant coefficients. The solution is approximated as

The operator $\mathbf {Q}(z)$ computes the scattering induced by the RI difference in the object slice at $z$ and is computed directly in the real space as

2.2 SSNP-based sequential IDT forward model

The sequential IDT, such as annular IDT [7], uses a single LED to illuminate the sample in each measurement. To derive the SSNP-based forward model, we first compute the illumination field and its axial derivative at $z=z_0$ as the initial condition to the SSNP equations. Each illumination field is a plane wave $\varphi (\boldsymbol {r})=\varphi _0 \exp [j(k_x^{in}\,x + k_y^{in}\,y + k_z^{in}\,z)]$ where $\boldsymbol {k}^{in} = (k_x^{in}, k_y^{in}, k_z^{in})$ is the wave vector. By omitting the constant factors $\varphi _0$ and $\exp (j\,k_z^{in}\,z_0)$, we have

By propagating the illumination field through the sample with the SSNP equations, we obtain the exit field from the sample at plane $z_n$. To compute the field through the microscope reaching the camera $\boldsymbol {\Phi }_\mathrm {cam}$, we first back-propagate the exit field from $z_n$ to the front focal plane ($z_\mathrm {focal}$), and then filter it by the microscope’s pupil function:

To summarize, our SSNP-based sequential IDT (SSNP-sIDT) forward model can be written in the following compact form

2.3 SSNP-based multiplexed IDT forward model

Multiplexed IDT uses several LEDs to illuminate the sample in each measurement [10]. Since the light from different LEDs are incoherent with each other, we can model the scattering process independently and compute the measured intensity as the sum of the intensities from all the LEDs in each measurement. Therefore, the forward model of the SSNP-based multiplexed IDT (SSNP-mIDT) is similar to the sequential IDT but with an additional incoherent sum over Eq. (12):

2.4 Inverse problem for SSNP-sIDT

Next, we formulate the inverse problem for reconstructing the 3D RI distribution using multiple intensity measurements taken from $L$ different LEDs $I_\mathrm {meas}^l$ indexed by $l=1, 2, \dots, L$. We formulate the sequential IDT reconstruction as the following optimization problem

We calculate the gradient using the automatic differentiation technique. Since the SSNP forward model consists of a chain of simple operations, the gradient of the data-fidelity term can be computed using the chain rule to perform “error back-propagation”. We define the gradient $\mathcal {G}$ of a complex function $w(v)$ as $\mathcal {G}(w,v)=\overline {\partial w / \partial v+\partial \overline {w}/\partial v}$ according to the Wirtinger calculus [20], where the $\overline {\cdot }$ operator means the complex conjugate. The chain rule of the gradient is

With the gradient of the data-fidelity term, we perform the reconstruction using the fast iterative shrinkage-thresholding algorithm (FISTA) to solve Eq. (14). The 3D TV regularization is implemented using an efficient proximal operator based on [15]. We initialize $n(\boldsymbol {r})$ as the background RI, $n_0$, and then update the estimation iteratively.

2.5 Inverse problem for SSNP-mIDT

For multiplexed IDT data, we have intensity measurements from $L$ illumination patterns, with each pattern containing $M$ LEDs. Correspondingly, the reconstruction is achieved by solving the following optimization problem:

To derive the gradient for the multiplexed data-fidelity term, we discuss a single LED pattern and omit the pattern index in the subsequent derivation for brevity. The data-fidelity term for the $l$th measurement is $\mathcal {L}' = \left \lVert \sqrt {\sum _{m=1}^M I_\mathrm {out}^m}-\sqrt {I_\mathrm {meas}} \right \rVert _2^2$, and the local gradient is

Overall, our SSNP-IDT algorithm is computationally efficient and requires $\sim 2\times$ time and memory cost as the BPM algorithm since the SSNP requires computing both the field and its axial derivative. We implemented the algorithm using custom CUDA code on Python to leverage the massively parallel processing power of GPU [22]. To reconstruct a single volume containing 600 $\times$ 600 $\times$ 150 voxels, the algorithm typically requires 20 iterations to converge and takes $\sim {120}\;\textrm{s }$ on a PC equipped with a GeForce RTX2070 GPU. A unified efficient algorithm for both the sequential and multiplexed IDT and an evaluation of the computational performance are provided in Supplement 1.

3. Results

3.1 Simulation: forward model accuracy evaluation

A major benefit of the SSNP model compared with the BPM model is its higher accuracy when computing scattered field under high-NA illumination [16]. To evaluate the SSNP-IDT model’s accuracy under high-NA illumination, we simulated sequential IDT measurements of a bead using the BPM and SSNP models and compare them against the Mie scattering theory. In this simulation, the bead has an RI of 1.02 and is immersed in air ($n_0=1$). The diameter of the bead is $6\lambda = {3.09}\;\mathrm{\mu}\textrm{m}$ where $\lambda = {0.515}\;\mathrm{\mu}\textrm{m}$ is the wavelength. In Fig. 2, we show the on-axis and an off-axis 0.9 NA illumination case. The microscope has a 0.9 NA objective lens. The discretization steps are $(\Delta x, \Delta y, \Delta z)=(\lambda /4, \lambda /4, \lambda /8)$. The simulated intensities from the BPM and SSNP IDT models and their difference maps from the Mie theory are shown. For the on-axis illumination case, both the SSNP and BPM models agree with the Mie theory with negligible errors. For the off-axis high-NA illumination case, SSNP retains high accuracy while the BPM model suffers from more severe errors. This study confirms that our SSNP-IDT forward model is highly accurate in simulating IDT measurements under high-NA illuminations. As a result, we will use the SSNP model for simulating the measurements in the subsequent numerical studies due to its higher computational efficiency and flexibility compared to the Mie theory.

 

Fig. 2. Simulation using the SSNP and BPM IDT models for a $6\lambda$ diameter bead ($n=1.02$) in air ($n_0=1$) using (a) on-axis illumination (NA = 0), and (b) off-axis illumination (NA = 0.9). The bead center is placed at the origin at $z=0$, and the observation plane is placed at $z=125\lambda$. The upper row shows the normalized intensity, and the lower row shows the difference comparing against the Mie theory.

Download Full Size | PDF

3.2 Simulation: reconstruction algorithm evaluation

To assess our SSNP-IDT reconstruction algorithm in the non-paraxial regime, we simulated high-NA annular IDT measurements using the SSNP-sIDT forward model and then performed reconstruction using the algorithm in Section 2.4. We simulated a $6\lambda$-diameter sphere with voxel size $(\lambda /4, \lambda /4, \lambda /8)$. To study the performance under different scattering strengths, we used two RI contrast values, including $\Delta n$=0.01 and 0.05. We simulated 8 intensity images from illuminations evenly distributed on a 0.89 NA ring, collected by a 0.9 NA objective (Fig. 3(a)).

 

Fig. 3. (a) An example image from the SSNP simulated intensity measurements for low and high RI sphere samples. (b) Error maps of the RI reconstructions from the SSNP, original BPM, and modified BPM models. The X-Y and X-Z cross sections through the center of each sphere are shown. (c) The cutlines along X and Z of the reconstructed and ground-truth RI through the center of each sphere. (d) The relative MSE of all the reconstructions.

Download Full Size | PDF

We performed reconstruction using both the SSNP, the original BPM [15], and modified BPM model [16] (see details in Supplement 1) in Fig. 3(b) and 3(c). Then we use the relative mean squared error (MSE) $\mathscr {E}$ to quantify the overall reconstruction quality in Fig. 3(d), which is defined as $\mathscr {E}=\frac {{\left \lVert {n_{GT}-\hat {n}}\right \rVert }^2}{{\left \lVert {n_{GT}-n_0}\right \rVert }^2},$ where $n_{GT}$ and $\hat {n}$ are the ground-truth and reconstructed RI, respectively. The original BPM over-estimates the RI values in both the weakly and strongly scattering cases. This is because the BPM forward model underestimates the phase modulation for the high-NA components [16]. The modified BPM model introduces an cosine obliquity factor to compensate for the under-estimated phase modulation, which improves the model accuracy at high angles [16]. However, this approximation results in under-estimation in the reconstructed RI values. Figure 3(d) quantitatively shows that the SSNP model can recover the RI distribution with better accuracy compared to the BPM and modified BPM models in both the weakly and strongly scattering cases. By inspecting the cross sections and cutlines shown in Fig. 3(b) and 3(c), we concluded that the reconstruction in all three dimensions are limited by the finite Fourier coverage provided by the illumination and imaging optics. Both the lateral resolution and the amount of axial elongation due to the “missing-cone” problem are comparable regardless of the model used.

3.3 Experiments: sequential IDT

To test our SSNP-sIDT reconstruction algorithm, we apply it to the annular IDT data from our previous publication in [7]. Briefly, our annular IDT system consists of a Nikon E200 microscope equipped with a ring LED (Adafruit, 1586 NeoPixel Ring). We used a $40\times$/$0.65$ NA (Nikon, CFI Plan Achromat) objective. The ring LED has 24 LEDs and is 60 mm in diameter. It is centered at the optical axis and placed approximately 35 mm away from the sample, which sets the illumination angle to be $\sim {40}^{\circ}$ and matches with the objective NA. Each LED illumination is approximated as a plane wave with a wavelength $\lambda = {515}\;\textrm{nm}$. We use the same self-calibration method as [7] to get the accurate direction of the wave vector.

In our reconstruction algorithm, the physical quantities on the XY plane, including the field, the field’s axial derivative, and the slices of the predicted RI values, are discretized on 2D grids. The grid spacing is $\Delta x=\Delta y= {0.1625}\;\mathrm{\mu}\textrm{m}$, which is chosen to be smaller than but close to $\frac {\lambda }{4 \mathrm {NA}}$ to ensure both accuracy and computational efficiency of the algorithm, and is the same as the effective pixel size on the camera.

First, we verify the reconstruction algorithm on unstained buccal epithelial cells immersed in purified water on a glass slide and covered by a coverslip. The dataset contains $24$ intensity images. The reconstructed volume is between -9.7 $\mu$m to 9.7 $\mu$m around the focal plane, and the FOV is 162.5 $\mu$m $\times$ 162.5 $\mu$m. We discretized the sample volume to $150$ slices, each with $1000 \times 1000$ pixels, with voxel size 0.1625 $\mu$m $\times$ 0.1625 $\mu$m $\times$ 0.129 $\mu$m. To perform the reconstruction with memory-limited hardware, we cropped the full FOV into four 576 $\times$ 576-pixel subregions, reconstructed them separately, and then stitched them together. Figure 4(a) shows the 3D rendering of the reconstructed RI distribution of the entire cell cluster volume. The 3D reconstruction allows easily discriminating cells at different depths. As shown in Fig. 4(b)-(c), our SSNP algorithm successfully reconstructs high-resolution features, such as cell boundaries, membrane, and native bacteria around the cells.

 

Fig. 4. Reconstruction of buccal epithelial cells. (a) 3D rendering of the full-FOV reconstruction. (b), (c) Zoomed-in views of XY (Top), XZ (Bottom) cross sections. The red ellipses show the native bacteria crowded regions. The yellow arrows highlight the sharp cell boundary features.

Download Full Size | PDF

Next, we apply our algorithm to a live C. elegans worm sample. Young adult C. elegans were mounted on 3% agarose pads in a drop of nematode growth medium (NGM) buffer. Glass coverslips were then gently placed on top of the pads and sealed with a 1:1 mixture of paraffin and petroleum jelly. The time series IDT measurements were taken at $\sim {10}\;\textrm{Hz}$ with 8 intensity images for each set. The reconstruction is performed in a volume of 97.5 $\mu$m $\times$ 390 $\mu$m $\times$ 38.6 $\mu$m centered at the focal plane. The volume is discretized to 600 $\times$ 2400 $\times$ 300 voxels with voxel size 0.1625 $\mu$m $\times$ 0.1625 $\mu$m $\times$ 0.129 $\mu$m. The reconstruction is more challenging since the worm contains complex and strongly scattering features. Figure 5(a) shows four consecutive volumes reconstructed by our SSNP algorithm, displayed as color-coded depth projections. Figure 5(b)-(c) shows two zoomed-in XY and YZ cross-sectional views around the buccal cavity and terminal pharyngeal bulb regions. Comparing with the SSNP model, the first-Born model underestimates the RI value and misses the low-frequency information. This effect is also reported in other multiple-scattering models for IDT reconstruction [8,14]. Additional zoom-in regions at different depths and time points are shown in Fig. 5(d)-(e), containing features including the anterior bulb, vulva, lipid droplets, and body wall muscle cells inside the worm sample.

 

Fig. 5. Reconstruction of a live C. elegans worm. (a) Color-coded depth projections of the whole worm reconstruction using SSNP model at different times. (b-c) XY and ZY cross-sections of the buccal cavity (top), part of the isthmus and terminal bulb (bottom). The reconstructions are using SSNP model and first-Born based model, respectively. (d) XY cross-sections of SSNP-IDT reconstruction at different time points of the same subregion and depth. (e) XY cross-sections of SSNP-IDT reconstruction at different depths and time points of the same subregion.

Download Full Size | PDF

3.4 Experiment: multiplexed IDT

Next, we test our SSNP-mIDT reconstruction algorithm to process data from our previous publication in [10]. Briefly, the multiplexed IDT system consists of a Nikon TE 2000U microscope equipped with a custom LED matrix. Each LED approximately generated a plane wave with a central wavelength $\lambda = {632}\;\textrm{nm}$. The measurements were taken with a 40$\times$ / 0.65 NA objective (Nikon, CFI Plan Achromat) and an sCMOS camera (PCO.Edge 5.5).

We apply our algorithm to a highly scattering sample of live C. elegans embryos. The multiplexed measurement used $16$ disjoint illumination patterns, each containing 6 LEDs. In total 96 LEDs are used in this experiment, corresponding to illumination NA ranging from $0.3$ to $0.575$. We reconstructed a $81.3 \;\mu\mathrm{m} \times 81.3\;\mu\mathrm{m} \times 19.0\;\mu\mathrm{m}$ volume, which was discretized to $500\;\mu\mathrm{m} \times 500\;\mu\mathrm{m} \times 120\;\mu\mathrm{m}$ voxels with voxel size $0.1625\;\mu\mathrm{m} \times 0.1625\;\mu\mathrm{m} \times 0.158\;\mu\mathrm{m}$. In Fig. 6, the embryo on the left is in the late three-fold (quickening) stage, and the one on the right is in the comma stage. From the color-coded view in Fig. 6(a), how the worms are folded can be clearly observed. From the single-depth cross-sections in Fig. 6(b), the morphological details of the cells’ outline, the buccal cavity, and the tail of the worm are reconstructed. Similar to aIDT, the SSNP model reconstructs more low-frequency phase information and higher RI contrast than the first-Born model. Besides, the TV regularization used in the SSNP model effectively suppresses the noise, which helps to provide a clean background and sharp embryo boundary comparing with the Tikhonov regularization used in the first-Born model. The high acquisition rate of multiplexed IDT enabled reconstructing the developing process and the moving behavior of the embryos, as shown in the time-series reconstructions. As compared to our previous single-scattering based linear multiplexed IDT reconstruction algorithm [10], SSNP-mIDT provides more robust demultiplexing capability that removes the fuzzy diffraction artifacts inside the embryos. Since the algorithm allows independent update from individual LEDs in each pattern, it greatly suppresses the cross-talk artifacts suffered by the linear model.

 

Fig. 6. Reconstruction of live C. elegans embryos. (a) Color-coded depth projection of SSNP-mIDT reconstructon at $t= {7.5}\;\textrm{s}, 15\;\textrm{s}, 30\;\textrm{s}$ (b) XY zoom-in cross sections at various times and depths using SSNP and first-Born models, which shows developing cells, buccal cavity and tail region of the embryos and worms, respectively.

Download Full Size | PDF

4. Conclusion and discussion

In conclusion, we developed a non-paraxial multiple-scattering model for 3D RI reconstruction from both sequential and multiplexed intensity diffraction tomography measurements. The SSNP model ensures high accuracy when imaging strongly scattering samples under high-angle illumination. We demonstrated the robustness of the model on weakly scattering cells and multiple-scattering dynamic C. elegans worm and embryos. With a unified reconstruction algorithm, our model is flexible for both the sequential and multiplexed IDT setups. This treatment for multiplexed measurement can also be applied to interferometric setups and speed up the acquisition for the ODT systems [1,16].

A major limitation of our current IDT systems is the limited angular coverage (up to 0.65 NA) that creates a significant amount of missing-cone artifacts and fundamentally limits the reconstruction quality, as shown by our numerical studies and experimental results. While several recent works alleviate the missing-cone artifacts by using higher NA (>1) optics and more than 10$\times$ more measurements [8,14], it sacrifices the FOV and acquisition speed, which prevented them from capturing the whole worm without mechanical FOV scanning and imaging dynamic biological processes. It may be possible to integrate the annular and multiplexed illumination strategies in these systems to speed up the acquisition while achieving better Fourier coverage.

It has been shown that 3D Fourier Ptychography models based on the first-Born approximation [23] and BPM [6] can further handle dark-field measurements, which allows considering high spatial frequency information beyond the objective NA in the 3D reconstruction. Since our SSNP model can be treated as an enhancement of the first-Born and BPM models, it is possible to extend our work to incorporate additional dark-field information for high-resolution IDT reconstructions, which will be considered in our future work. Alternatively, several advanced reconstruction frameworks for IDT have been recently developed and shown promising results to greatly suppress the missing-cone artifacts, including end-to-end supervised learning [24], untrained neural network [9], model-based reconstruction with a deep denoiser [19], neural fields based reconstruction [25], and physics-informed neural network [26]. Notably, our SSNP forward model has been previously applied to generate a large-scale training data set to train a supervised learning model [24]. It is further possible to directly incorporate our SSNP model into the model-based learning strategies [9,19,25] in the future.