Highfidelity optical diffraction tomography of multiple scattering samples  Light: Science & Applications

Abstract.We propose an iterative reconstruction scheme for optical diffraction tomography that exploits the splitstep nonparaxial (SSNP) method as the forward model in a learning tomography scheme. Compared with the beam propagation method (BPM) previously used in learning tomography (LTBPM), the improved accuracy of SSNP maximizes the information retrieved from measurements, relying less on prior assumptions about the sample. A rigorous evaluation of learning tomography based on SSNP (LTSSNP) using both synthetic and experimental measurements confirms its superior performance compared with that of the LTBPM. Benefiting from the accuracy of SSNP, LTSSNP can clearly resolve structures that are highly distorted in the LTBPM. A serious limitation for quantifying the reconstruction accuracy for biological samples is that the ground truth is unknown. To overcome this limitation, we describe a novel method that allows us to compare the performances of different reconstruction schemes by using the discrete dipole approximation to generate synthetic measurements. Finally, we explore the capacity of learning approaches to enable data compression by reducing the number of scanning angles, which is of particular interest in minimizing the measurement time. Introduction.Quantitativephase imaging (QPI) enables the measurement of the phasecontrast information of transparent samples such as biological cells. QPI contrast is generated from the refractive index (RI) contrasts within and around a sample. Because this contrast mechanism is endogenous, quantitativephase information does not require external labeling, such as immunostaining, which may perturb the sample. QPI contains the coupled information of sample thickness and RI contrast. Optical diffraction tomography (ODT) provides the 3D RI distribution of a sample by combining multiple 2D QPI measurements from various illumination angles1,2. Reconstructed tomograms provide structural information that has been extensively utilized to study hematology3,4, morphological parameters5, and biochemical information6 which are summarized in several review papers2,7,8,9. In ODT, the way in which multiple 2D measurements are combined into unified 3D information is critical. Under the assumption of a weakly scattering sample, the Wolf transform10 has been widely used. Depending on how the 2D projections are processed, we obtain either the Born or Rytov approximations for the Wolf transform1. Each method has its limitations11, but the Rytov approximation is known to be more appropriate than the Born approximation for many biological applications12. However, when a sample is thicker and more complex, the Rytov approximation is no longer valid. This limits the usefulness of ODT for imaging complex samples.Recently, methods have emerged to overcome the limitations of the Born and Rytov approximations by taking multiple scattering into account13,14,15,16,17,18,19,20. It was shown using Mie theory20 that learning tomography (LT)14,21, an approach that exploits the beam propagation method (BPM) as the forward model to capture multiple scattering, has superior performance compared with that of the conventional imaging method based on the Rytov approximation. We refer to it as LTBPM. LT uses the forward model of dividing 3D samples into multislices followed by slicebyslice propagations. Due to the multislice modeling of forward models by LT, the resulting structure is similar to that of a neural network, and we can use the error backpropagation algorithm to calculate the gradient. The BPM consists of two steps: nonparaxial diffraction followed by phase modulation. The diffraction step used in the BPM assumes \(k_0n\left( {x,y,z} \right) \approx k_0n_0\) , where k0 is the freespace wavenumber, n0 is the RI of the medium, and n(x,y,z) represents RI variations. In addition, the phase modulation steps use a distance, dz/cosθ, to modulate the phase throughout propagation, given the propagation step (dz) and the illumination angle (θ). However, for thicker and more complex samples, as light propagates through the samples, multiple diffracted beams of light are generated, and it is not valid to use one single value, dz/cosθ, to represent optical path lengths. This deviation from the fixed distance, dz/cosθ, increases with increasing the illumination angle due to the nature of the cosine function22.In this paper, we show that the accuracy of LT reconstructions of a 3D object is increased when we use the splitstep nonparaxial (SSNP) method rather than the BPM. We refer to it as LTSSNP. The SSNP method exploits not only the field but also the derivative of the field along the optical axis to model the propagation23,24. While the BPM requires this approximation,\(k_0n\left( {x,y,z} \right) \approx k_0n_0\) , to decouple diffraction from phase modulation, SSNP does not require the approximation, benefiting from propagating the derivative of the field at the same time. Phase modulation affects the derivative and is used concurrently in the next step of the diffraction calculation. LTSSNP uses the same iterative scheme used in LTBPM. To fairly assess LTSSNP and compare it with the LTBPM, synthetic measurements are generated using Mie theory and the discrete dipole approximation (DDA). For spherical and cylindrical objects, Mie theory provides the analytical solution to the Helmholtz equation25. Therefore, the solution of Mie theory takes into account multiple scattering. Here, we also use the DDA to simulate light scattering by an arbitrarily shaped sample to generate more complex synthetic data. The DDA is a general method for calculating the scattering and absorption caused by an arbitrarily shaped sample represented by finite discrete dipoles26. These dipoles react not only to incident light but also to one another, which places the resulting fields under high orders of scattering. It has been shown that the DDA works well for samples whose RI values fairly match those of the surroundings, such as biological cells in a liquid medium27. Therefore, we use Mie theory for multiple cylinders and the DDA for a cell phantom, as well as a cluster of 15 red blood cells (RBCs). After generating synthetic measurements by using either Mie theory or the DDA, the LTBPM and LTSSNP are used to reconstruct the 3D RI of each sample, and the accuracy of each reconstruction is evaluated quantitatively.In this analysis, we include an investigation of the performance of each algorithm with respect to regularization. The iterative reconstruction scheme used for both the LTBPM and LTSSNP minimizes a cost function that comprises two terms: data fidelity and regularization. The data fidelity term is defined by whether the forward model applies either the BPM or the SSNP, and the regularization term introduces prior knowledge about the sample characteristics such as edge sparsity and nonnegativity. The relative importance of the two terms in the cost function is controlled by the regularization parameter. We compare the LTBPM and LTSSNP by using varying regularization parameters with the goal of minimizing the influence of the regularization term so that the results are primarily based on the forward model rather than on prior knowledge. For the simulations described, we confirm that LTSSNP shows lower dependency on the regularization parameter due to the accuracy of SSNP. In other words, the use of a more accurate forward model permits LTSSNP to extract more information from the measurements and to rely less on regularization. More importantly, for highly aggregated samples subject to significant multiple scattering, LTSSNP allows individual objects and structures to be clearly distinguished, while this observation cannot be made when using the LTBPM.We validate the proposed method by using experimental ODT data from a yeast cell and from HCT116 human colon cancer cells. To image biological cells with fine details, it is critical to reduce the influence of the regularization term, as high regularization not only smooths out the imaging artifacts but also useful information, leading to deterioration in the quality of the reconstruction. Tomograms of a yeast cell reconstructed by using LTSSNP show successful results with high quality even with a very low regularization parameter, while the LTBPM fails to recover fine details within and around the cells. In the case of experimental measurements of biological cells, the true RI distribution is not known, which prevents the direct assessment of the accuracy of the various ODT methods. To overcome this issue, we generate two sets of semisynthetic measurements by using the DDA for each of the RI reconstructions from the LTBPM and LTSSNP. A comparison of the discrepancies between the semisynthetic and experimental measurements reflects the proximity of each solution to the real RI values.Finally, we explore the capacity of LTSSNP to produce accurate reconstructions with a reduced number of illumination angles28,29. This is of particular interest because the number of scanning angles is directly related to the measurement time. A comparison of each reconstruction method for a varying number of scanning angles indicates that learning approaches provide a dramatic improvement over conventional methods. Overall, the more accurate forward model used in LTSSNP translates to excellent results even with low regularization and a small number of illumination angles.Results.In this section, we compare the LTBPM and LTSSNP, which belong to the same family of LT reconstruction schemes, except for the forward models, namely, the BPM and the SSNP, respectively. LT minimizes the cost function, which consists of two terms as follows:$$\widehat {\bf{x}} = \arg \min _{{\bf{x}} \in P}\frac{1}{{2L}}\mathop {\sum}\limits_{l = 1}^{L} {\mathbf{y}}_{K}^{(l)}  {\mathbf{S}}_{K}^{(l)}({\bf{x}})_{2}^{2} + \tau R({\bf{x}})$$ (1) where the first term is the data fidelity term and R is the 3D total variation (TV)30 regularization term to impose edge sparsity on the solution. The relative importance between two terms is controlled by the regularization parameter, τ. \({\mathbf{y}}_K^{\left( l \right)} \in {\Bbb C}^M\) denotes the experimental measurements at the Kth slice for each illumination angle l, and L is the total number of angles. \({\mathbf{S}}_K^{\left( l \right)}\left( {\mathbf{x}} \right)\) represents the estimate by a forward model (either the BPM or the SSNP) at the Kth slice, which is the last slice of the volume, to be compared with \({\mathbf{y}}_K^{\left( l \right)}\) given a current solution, \({\mathbf{x}} \in {\Bbb R}^N\) . \(P \in {\Bbb R}^N\) is a convex set that imposes a nonnegativity constraint. In the supplementary section, we describe the calculation of the gradient for SSNP. Once we calculate the gradient of the data fidelity term in Eq. (1), the optimization scheme uses the fast iterative shrinkagethresholding algorithm (FISTA)31 as explained in ref. 21 for 3D isotropic TV regularization, with eight randomly chosen angles in each iteration.Multiple cylinders by using Mie theory.We applied the LTBPM and LTSSNP on a highly scattering simulated sample consisting of a 3 × 3 grid of cylinders. Each cylinder is 6 μm in diameter with an RI of 1.05 immersed in air. The centertocenter distance is 9 μm. We varied the regularization parameter to investigate the accuracy of the forward model for each algorithm. The results are presented by mapping the difference between the reconstructed tomogram for each method and the known solution, as shown in Fig. 1a. The LTBPM shows many artifacts inside the cylinders and smearing of the RI in the region between the cylinders. These artifacts of the forward model cannot be eliminated even by increasing the regularization parameter. Regularization only smooths out the overall reconstruction. In contrast, LTSSNP clearly distinguishes each cylinder without interstitial artifacts even with the weakest regularization parameter tested, that is, 0.25τ = 0.01. Interestingly, increasing the regularization parameter to 4τ = 0.16 reduces the reconstruction quality when using the LTSSNP algorithm. The total Error, which is defined as follows:$$Error\,({\bf{x}}_{\mathrm{recon}},{\bf{x}}_{\mathrm{true}}) = \frac{{{\bf{x}}_{\mathrm{recon}}  {\bf{x}}_{\mathrm{true}}^2}}{{{\bf{x}}_{\mathrm{true}}^2}}$$ (2) was also calculated as a function of the iteration number, as shown in Fig. 1b. xrecon is the reconstructed RI contrast from the medium RI, and xtrue is the ground truth RI contrast. Figure 1b displays the Error plots of the LTBPM and LTSSNP by using the regularization parameter that produced the lowest Error value for each algorithm: 4τ = 0.16 for the LTBPM and τ = 0.04 for LTSSNP. This analysis quantitatively confirms the better accuracy of LTSSNP. In the case of multiple cylinders, it is critical to model distortions in the wavefront (phase modulation) introduced by the precedent samples, which determine the illumination on subsequent samples. We further analyzed this scenario by varying the number of layers in multiple cylinders and summarized the results in the supplementary material.Fig. 1: Reconstruction results of cylinders using the LTBPM and LTSSNP for various regularization parameters (τ = 0.04).a Difference maps between reconstructions from the LTBPM/LTSSNP and the ground truth (reconstruction—truth). b Plots of the Error of the LTBPM and LTSSNP by using the regularization parameter that produced the minimum Error valueFull size image .RBC cluster using discrete dipole approximation.To investigate the performance of each algorithm with highly scattering samples in 3D, we performed a similar test on a simulated cluster of RBCs. The shape of a single RBC is sketched in Fig. 2a, while the organization of the cluster is shown in Fig. 2b. Reconstructions were performed by using various regularization parameters; for each algorithm, we show only the reconstruction by using the regularization parameter that gives the lowest Error: 8τ = 0.2 for the LTBPM and τ = 0.025 for LTSSNP. In Figs. 2c, 3, different slices (xy, yz, and xz) of the 3D RI distributions resulting from each method are presented, along with the difference map with respect to the ground truth. Both the LTBPM and LTSSNP show better reconstructions compared with reconstructions based on the Rytov approximation, which is expected since Rytov does not consider multiple scattering. By comparing the LTBPM and LTSSNP, we can see that the RI tomogram resulting from LTSSNP shows clearer and more accurate reconstructions of each RBC, producing homogeneous RI distributions within each RBC.Fig. 2: Reconstruction results of a RBC cluster using the LTBPM and LTSSNP.a Parameters used to define the shape of the RBC. b 3D rendering of the RBC cluster, which consists of 15 identical RBCs. c 3D RI and difference maps for a cluster of 15 RBCs. Top to bottom: ground truth, Rytov approximation, LTBPM, and LTSSNP (τ = 0.025). Left: xy, yz, and xy slices of the 3D RI. Right: difference maps between the reconstructed RI and the ground truthFull size image Fig. 3: Reconstruction results of a cell phantom by using Rytov, the LTBPM, and LTSSNP at four different z planes (τ = 0.025).a Left column: the yz slice, with the zpositions of the yx slices indicated by using dashed lines. Second through fifth columns: yx slices at the positions indicated in the left column. b Magnification of the regions indicated by the boxes in part (a)Full size image .Cell phantom using discrete dipole approximation.To evaluate the LTBPM and LTSSNP algorithms on a sample whose RI values are not homogeneous and which contains fine details, we generated a synthetic cell phantom. The phantom contains four different RI values corresponding to the cytoplasm, nucleus, nucleolus, and lipids32. Synthetic measurements were made by using the DDA in the same manner as for the RBCs. Again, we present for each algorithm only the results obtained by using the regularization parameter that produced the lowest Error value (4τ = 0.1 for the LTBPM and τ = 0.025 for LTSSNP). The reconstruction results are shown in Fig. 3a. Figure 3b displays the magnified regions, which are indicated by the yellow boxes in Fig. 3a. The reconstructions obtained with the LTBPM show that the cytoplasmic regions are highly distorted, similar to the artifacts observed inside the cylinders and RBCs. More importantly, the small lipids indicated by the red arrows are hardly distinguishable due to the inaccuracy of the BPM. By contrast, the LTSSNP not only distinguishes the shapes of fine structures but also correctly positions them along the optical axis.Experimental validation using a yeast cell.To validate the relative performances of the LTBPM and LTSSNP on experimental data, we acquired ODT images of a yeast cell. Again, we evaluated different regularization parameters for the reconstructions obtained by using the LTBPM and LTSSNP. Figure 4 shows the reconstruction results for a slice close to the image plane. For both the LTBPM and LTSSNP, the high regularization parameter 4τ = 0.1 results in too much smoothing, and it becomes difficult to resolve fine details. Therefore, it is necessary to reduce the regularization parameter. However, in the case of the LTBPM, lowering the value of the regularization parameter introduces artifacts similar to those present in the simulation results in the previous section. By contrast, the LTSSNP can reconstruct fine details without introducing strong artifacts. Therefore, we used τ = 0.025 for the LTBPM and τ/4 = 0.00625 for LTSSNP and further analyzed the sample for different z planes, as shown in Fig. 5. Since we used higher regularization for the LTBPM, we can clearly see that images tend to be smoothed out and fine details are lost, as indicated by the red arrows in Fig. 5. By contrast, the LTSSNP reveals structures that are not observable in the Rytov and LTBPM reconstructions. In addition, even with the higher regularization, the LTBPM still shows several artifacts, as indicated by the black arrows in Fig. 5.Fig. 4Reconstruction results of a yeast cell by using the LTBPM and LTSSNP for various regularization parameters (τ = 0.025)Full size image Fig. 5Reconstruction results of a yeast cell by using Rytov, the LTBPM, and LTSSNP at three different z planes (τ = 0.025)Full size image A serious limitation for quantification of the reconstruction accuracy for real biological samples such as this yeast cell is the fact that the true RI distribution is unknown. However, we were able to further evaluate the differences between the LTBPM and the LTSSNP by using the semisynthetic measurements generated by using the DDA. While we generated synthetic measurements by using synthetic samples in all previous cases, the RI reconstructions obtained by using the LTBPM and the LTSSNP served as samples for the DDA to generate semisynthetic measurements in this case. The projection error—the difference in phase information between the experimental data and these simulated measurements—reflects how close the solution is to the true RI distributions, as shown in Fig. 6a. Figure 6b maps the 2D projection error for two randomly selected angles as well as the average across the full set of angles for each algorithm. In the case of the LTBPM, differences are clearly observed when compared with the LTSSNP, which shows remarkable consistency with the experimental measurements. We quantified the mean projection error (radians/pixel) for each and used this metric to quantify the accuracy of the LTSSNP compared with that of the LTBPM (Fig. 6c). The average projection error across all angles was 65% lower for the LTSSNP than for the LTBPM.Fig. 6: Evaluation of LT algorithms by using semisynthetic error estimation.a Overall scheme of semisynthetic measurement generation by using the DDA. b Phasedifference maps for two randomly selected angles and the average for all angles. The color bar is in radians. c Calculation of the projection error in retrievedphase information from experimental measurements and semisynthetic dataFull size image .Data compression demonstrated on experimental data—HCT116 cells.The tomographic reconstruction based on the Wolf transform and the Rytov approximation directly maps multiple 2D measurements into the 3D Fourier space. Therefore, any missing information in measurements directly deteriorates the final reconstruction. However, the LTSSNP is an iterative reconstruction scheme. The iterative reconstruction begins with an initial guess (usually based on the Rytov approximation), and the initial solution is updated based on the calculated error gradient by using the forward model. In addition, prior knowledge about the sample is imposed on the current guess during the iterative process. Therefore, even if the measurements are underdetermined due to missing measurements, the learning approaches can fill in some of the missing information. This idea was validated by reducing the number of illumination angles used for each method. The experimental data used for this investigation were ODT images of a pair of HCT116 human colon cancer cells. These cancerous epithelial cells contain information in small structures relative to the size of the cell and highlight the importance of reconstructions that can capture these fine details. Reconstructions were performed by using Rytov, linear tomography20, and the LTSSNP by using different numbers of projection angles (45, 24, 12, and 4) uniformly spaced in the range from 0 to 360°. The linear tomography method uses the same iterative reconstruction scheme as the LTSSNP, except with single scattering as the forward model. For the quantitative analysis, we also compare the structural similarity index (SSIM)33 for reconstructions from compressed measurements with the full measurement case, namely, 360 angles at the focal plane. The results, plotted in Fig. 7, show a dramatic improvement in the reconstruction quality for linear tomography and the LTSSNP because the two methods iteratively fill up empty components introduced from missing measurements but using different forward models. In the case of the HCT116 cells, Rytov produces fairly good reconstructions that reveal intracellular structures with 360 full projections, despite the underestimation due to the missingcone problem. The Rytov reconstructions, on the other hand, rapidly deteriorate as the number of illumination angles decreases. Compared with Rytov and linear tomography, the LTSSNP is more robust in the number of projections, providing reconstructions with only four scanning angles with nearly the same quality as reconstructions by using the full 360angle data, as confirmed by the SSIM in Fig. 7b. We believe that the LTSSNP can benefit from both the iterative scheme and an accurate forward model. In addition, we further tested the compression using the cell phantom, which has higher RI contrasts; the results have been added to the supplementary material.Fig. 7: Data compression demonstrated on HCT116 cells.a Reconstruction results of HCT116 human colon cancer cells by using Rytov, linear tomography, and the LTSSNP for a downsampled number of scanning angles. b Structural similarity index plot with respect to the reconstruction by using each method with full measurementsFull size image .Discussion.In this study, we have proposed a new tomographic reconstruction algorithm, the LTSSNP, which is based on the SSNP forward model, for imaging complex highly scattering samples with fine details. By benefiting from the accuracy of the SSNP, the LTSSNP extracts a maximum amount of information from measurements rather than relying on prior assumptions and generalizations about the sample structure. The LTSSNP was quantitatively evaluated and compared with the previous algorithm, the LTBMP, by using synthetic measurements. These synthetic measurements with a known solution were generated by using Mie theory for multiple cylinders, and the DDA for an arbitrarily shaped cluster of RBCs and a cell phantom.In the case of multiple cylinders, the LTSSNP shows clear reconstruction of each sample without introducing artifacts. The more interesting point is that the LTSSNP does not require strong regularization. This is because the SSNP forward model is accurate enough that regularization is not necessary to compensate for poor data fidelity, while the LTBPM could not properly carry out the reconstruction even with high regularization. For the RBC cluster in 3D, the LTSSNP returns more homogeneous distributions even with a lower value of the regularization parameter than that of the LTBPM. This fact is critical when imaging complex samples because too much regularization smooths out fine structures and makes them impossible to resolve. The cell phantom simulation confirms the performance of the LTSSNP on a sample with highresolution information. The LTSSNP is more accurate and permits the use of a lower regularization parameter, which allows details of the 3D refractive index to be identified without artificially being smoothed out by regularization.Importantly, the added capabilities of the LTSSNP are dramatic for imaging biological samples containing information across many scales, as confirmed by applying it to tomographic images of a yeast cell. The reconstructed tomograms by using the LTSSNP clearly reveal structures that are not observable in the case of Rytov and the LTBPM. Semisynthetic measurements based on the RI reconstructions of the LTBPM and the LTSSNP numerically validate the accuracy of the LTSSNP reconstructions. The averaged phasedifference map represents how close the reconstruction using each method is to the real sample. In contrast to the averaged phasedifference map of the LTBPM, which produces many discrepancies inside the sample, that of the LTSSNP shows consistency with the experimental measurements. The numerical evaluation shows that the LTSSNP produces a 65% reduction in the projection error compared with that of the LTBPM.Furthermore, we explored the capacity of learning approaches to enable data compression by reducing the number of scanning angles. The LTSSNP shows a dramatic improvement in image quality by using a small number of illumination angles when compared with the conventional direct inverse method by using the Rytov approximation. Even with a low number of projections, the LTSSNP benefits from its weak dependency on the regularization parameter.Materials and methods.Simulation.We used Mie theory to derive the field scattered by multiple cylinders (2D)34. A total of 101 illumination angles uniformly distributed between −45° and 45° were used. To perform a deeper assessment of the LTBPM and LTSSNP algorithms, we also tested on synthetic measurements in arbitraryshaped samples: an RBC cluster and a cell phantom.For RBC simulations, the discrete dipole approximation26,35 was applied to an RBC cluster, in which the surface of each RBC is defined by using the following equation:$$\rho ^4 + 2S\rho ^2z^2 + z^4 + P\rho ^2 + Qz^2 + R = 0$$ (3) where ρ is the radius in cylinder coordinates (ρ2 = x2 + y2) and S, P, Q, and R are parameters derived from d, h, b, and c shown in Fig. 2a, respectively. In this paper, d, h/d, b/d, and c/d were set to 7.7 μm, 0.3542, 0.1752, and 0.6196, respectively, as suggested in ref. 36. We refer interested readers to previous studies36,37 for a more complete presentation of the DDA simulation of an RBC. By using a single simulated RBC, a cluster consisting of 15 identical RBCs was generated, as shown in Fig. 2b. In addition, we generated a synthetic cell phantom with four different RI values corresponding to the cytoplasm, nucleus, nucleolus, and lipids32. To derive the scattered field from the cluster and the cell phantom, samples were scanned by using 40 uniformly distributed illumination angles on a circle with an incident angle of 45°. For every simulation mentioned above, the sample with an RI of n was immersed in air, and the wavelength used was 600 nm. This is equivalent to a case in which the RI of the medium is n0 and the sample with an RI of n × n0 is illuminated at a wavelength of 600 × n0 nm. The number of dipoles per wavelength for both simulations was set to 12. Table 1 summarizes the numerical and experimental parameters used for the simulations as well as for the experiments.Table 1 Reconstruction parametersFull size table .Experiments.The experiments were performed by using a conventional optical diffraction tomography configuration in which a spatial light modulator was used to control the illumination angle. A total of 360 holograms were recorded for each sample in a circular pattern with 1° resolution at an incidence angle of 35°. Additional details about the optical setup and sample preparation are provided in the supplementary section.Semisynthetic simulation.The semisynthetic measurements were calculated by using the reconstruction results acquired from the LTBMP and the LTSSNP as samples for the DDA. The size of the dipole was set to \(\frac{\lambda }{{12{{n}}_0}} = 0.033\,{\mathrm{nm}}\) , where λ = 0.532 nm is the wavelength of the laser and n0 = 1.338. Both the values were set from the values used in the experiments. The grid size of the reconstructions from the LTBPM and the LTSSNP was 99 nm. The reconstruction results were interpolated to a grid, one pixel of which was the size of a dipole. Then, we quantized the RI values by using the following equation: \({\mathrm{round}}(\frac{{n_{{\mathrm{recon}}}}}{{n_0}}\times1000)/1000\) , where nrecon denotes the reconstructed RI values. Simulations were performed for 160 nonoverlapping angles, which were calculated from the experiments.Reconstruction algorithm.We implemented the algorithms by using custom scripts in MATLAB R2018a (MathWorks Inc., Natick, MA, USA) on a desktop computer (Intel Core i76700 CPU, 3.4 GHz, 32 GB of RAM). To accelerate the computation, a graphicprocessing unit (GPU, GeForce GTX 1070) with custommade functions based on the compute unified device architecture (CUDA) was utilized. The gradient, calculated from a data fidelity term, D(\(\bf{x}\) ), was ∂D(\(\bf{x}\) )/∂\(\bf{x}\) , the amplitude of which is proportional to the amplitude of D(\(\bf{x}\) ). The LTBPM and the LTSSNP use different data fidelity terms. The LTBPM calculates the difference in the fields u(x,y,z). In contrast, the LTSSNP requires differences in both u(x,y,z) and its derivative du(x,y,z)/dz. Therefore, calibration of the optimization parameters between the methods is necessary to make the LTBPM and the LTSSNP use similar optimization parameters. The FISTA requires two parameters: step size (γ) and regularization parameter (τ). The calibration of those parameters can be performed by calculating the ratio C between \(u\left( {x,y,z} \right)_2^2\) and \(u\left( {x,y,z} \right) + du\left( {x,y,z} \right)/dz_2^2\) . We approximated this as the average value of (1+1ikz)2 for the illumination kzs, which corresponds to a case in which u(x,y,z) is replaced with a planar wave, \(e^{1i(k_xx + k_yy + k_zz)}\) . Therefore, the LTBPM, which uses the parameters γ and τ, can be directly compared with the LTSSNP, which uses the parameters γ/C and τ × C. For convenience, we labeled the figures according to the parameters used for the LTBPM. The actual parameter values for the LTSSNP can be easily calculated given C, which is provided in Table 1. The total number of iterations used in the FISTA is also provided in Table 1. Twenty iterations were used for the TV optimization step in all cases.Overall scheme of the learning tomography.Both algorithms (LTBPM and LTSSNP) start from measured electric fields (including both amplitude and phase information) from the holographic data. An initial guess of the RI distributions is obtained by using the Rytov tomographic reconstruction method. By using either the BPM or the SSNP as the forward model, the scattered field is estimated given the planewave illumination propagating through this initial guess. The square of the difference between the estimated and the measured fields is the cost function, which is minimized by adjusting the index values contained in the forward model through the FISTA. At the same time, an intermediate step of regularizations such as smoothness and nonnegativity is included. This process is repeated until the total cost function converges.Splitstep nonparaxial method.In this section, we briefly describe the SSNP in 3D23,24, which is the physical forward model used in the LTSSNP. Bhattacharya and Sharma38 implemented this method by using a matrix formalism for wave propagation in 3D. Here, we describe a fast Fourier transform implementation for more efficient use of memory.The propagation of a scalar wave u(x,y,z) through a medium n(x,y,z) in 3D must satisfy the following wave equation:$$\left( {\frac{{\partial ^2}}{{\partial x^2}} + \frac{{\partial ^2}}{{\partial y^2}} + \frac{{\partial ^2}}{{\partial z^2}}} \right)u\left( {x,y,z} \right) + k_0^2n^2\left( {x,y,z} \right)u\left( {x,y,z} \right) = 0$$ (4) where k0 = 2π/λ is the freespace wavenumber for a given wavelength λ in a vacuum. Eq. (4) can be written in matrix form$$\frac{\mathrm{d}{{\mathbf{v}}(x,y,z)}}{\mathrm{dz}} = {\mathbf{H}}(x,y,z){\mathbf{v}}(x,y,z)$$ (5) where,$${\mathbf{v}}(x,y,z) = \left[ {\begin{array}{{20}{c}} {u(x,y,z)} \\ {\frac{{u(x,y,z)}}{{\partial z}}} \end{array}} \right]$$ (6) and$${\mathbf{H}}(x,y,z) = \left[ {\begin{array}{{20}{c}} 0 & 1 \\ {  \left( {\frac{{\partial ^2}}{{\partial x^2}} + \frac{{\partial ^2}}{{\partial y^2}} + k_0^2n^2(x,y,z)} \right)} & 0 \end{array}} \right]$$ (7) When we consider an inhomogeneous sample immersed in a homogeneous medium, n0, it is possible to split the matrix H into two terms that correspond to diffraction and phase modulation. Note that no approximation is assumed up to this point. We refer interested readers to the supplementary section for a detailed explanation. References.1. Born, M. & Wolf, E. Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. (Cambridge University Press, 1999).2. Lee, K. et al. Quantitative phase imaging techniques for the study of cell pathophysiology: from principles to applications. Sensors 13, 4170–4191 (2013).Article. Google Scholar.3. Yoon, J. et al. Labelfree characterization of white blood cells by measuring 3D refractive index maps. Biomed. Opt. Express 6, 3865–3875 (2015).Article. Google Scholar.4. Park, Y. et al. Refractive index maps and membrane dynamics of human red blood cells parasitized by Plasmodium falciparum. Proc. Natl Acad. Sci. USA 105, 13730–13735 (2008).ADS.Article. Google Scholar.5. Kim, K. et al. Threedimensional labelfree imaging and quantification of lipid droplets in live hepatocytes. Sci. Rep. 6, 36815 (2016).ADS.Article. Google Scholar.6. Cooper, K. L. et al. Multiple phases of chondrocyte enlargement underlie differences in skeletal proportions. Nature 495, 375–378 (2013).ADS.Article. Google Scholar.7. Kim, K. et al. Optical diffraction tomography techniques for the study of cell pathophysiology. J. Biomed. Photonics Eng. 2, 020201 (2016). Google Scholar.8. Jin, D. et al. Tomographic phase microscopy: principles and applications in bioimaging [Invited]. J. Opt. Soc. Am. B 34, B64–B77 (2017).Article. Google Scholar.9. Park, Y., Depeursinge, C. & Popescu, G. Quantitative phase imaging in biomedicine. Nat. Photonics 12, 578–589 (2018).ADS.Article. Google Scholar.10. Wolf, E. Threedimensional structure determination of semitransparent objects from holographic data. Opt. Commun. 1, 153–156 (1969).ADS.Article. Google Scholar.11. Slaney, M., Kak, A. C. & Larsen, L. E. Limitations of imaging with firstorder diffraction tomography. IEEE Trans. Microw. Theory Tech. 32, 860–874 (1984).ADS.Article. Google Scholar.12. Sung, Y. et al. Optical diffraction tomography for high resolution live cell imaging. Opt. Express 17, 266–277 (2009).ADS.Article. Google Scholar.13. Tian, L. & Waller, L. 3D intensity and phase imaging from light field measurements in an LED array microscope. Optica 2, 104–111 (2015).Article. Google Scholar.14. Kamilov, U. S. et al. Learning approach to optical tomography. Optica 2, 517–522 (2015).Article. Google Scholar.15. Kamilov, U. S. et al. A recursive Born approach to nonlinear inverse scattering. IEEE Signal Process. Lett. 23, 1052–1056 (2016).ADS.Article. Google Scholar.16. Liu, H. Y. et al. SEAGLE: Sparsitydriven image reconstruction under multiple scattering. IEEE Trans. Comput. Imaging 4, 73–86 (2018).MathSciNet.Article. Google Scholar.17. Soubies, E., Pham, T. A. & Unser, M. Efficient inversion of multiplescattering model for optical diffraction tomography. Opt. Express 25, 21786–21800 (2017).ADS.Article. Google Scholar.18. Lim, J. et al. Beyond BornRytov limit for superresolution optical diffraction tomography. Opt. Express 25, 30445–30458 (2017).ADS.Article. Google Scholar.19. Pham, T. A. et al. Versatile reconstruction framework for diffraction tomography with intensity measurements and multiple scattering. Opt. Express 26, 2749–2763 (2018).ADS.Article. Google Scholar.20. Lim, J. et al. Learning tomography assessed using mie theory. Phys. Rev. Appl. 9, 034027 (2018).ADS.Article. Google Scholar.21. Kamilov, U. S. et al. Optical tomographic image reconstruction based on beam propagation and sparse regularization. IEEE Trans. Comput. Imaging 2, 59–70 (2016).MathSciNet.Article. Google Scholar.22. Bao, Y. J. & Gaylord, T. K. Clarification and unification of the obliquity factor in diffraction and scattering theories: discussion. J. Opt. Soc. Am. A 34, 1738–1745 (2017).ADS.Article. Google Scholar.23. Sharma, A. & Agrawal, A. Nonparaxial splitstep finitedifference method for beam propagation. Opt. Quantum Electron. 38, 19–34 (2006).Article. Google Scholar.24. Sharma, A. & Agrawal, A. New method for nonparaxial beam propagation. J. Opt. Soc. Am. A 21, 1082–1087 (2004).ADS.Article. Google Scholar.25. Mie, G. Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen. Ann. der Phys. 330, 377–445 (1908).ADS.Article. Google Scholar.26. Yurkin, M. A. & Hoekstra, A. G. The discretedipoleapproximation code ADDA: capabilities and known limitations. J. Quant. Spectrosc. Radiat. Transf. 112, 2234–2247 (2011).ADS.Article. Google Scholar.27. Yurkin, M. A. et al. Systematic comparison of the discrete dipole approximation and the finite difference time domain method for large dielectric scatterers. Opt. Express 15, 17902–17911 (2007).ADS.Article. Google Scholar.28. Kaganovsky, Y. et al. Compressed sampling strategies for tomography. J. Opt. Soc. Am. A 31, 1369–1394 (2014).ADS.Article. Google Scholar.29. Brady, D. J. et al. Compressive tomography. Adv. Opt. Photonics 7, 756–813 (2015).ADS.Article. Google Scholar.30. Beck, A. & Teboulle, M. Fast gradientbased algorithms for constrained total variation image denoising and deblurring problems. IEEE Trans. Image Process. 18, 2419–2434 (2009).ADS.MathSciNet.Article. Google Scholar.31. Beck, A. & Teboulle, M. A fast iterative shrinkagethresholding algorithm for linear inverse problems. SIAM J. Imaging Sci. 2, 183–202 (2009).MathSciNet.Article. Google Scholar.32. Müller, P., Schürmann, M. & Guck, J. ODTbrain: a Python library for fullview, dense diffraction tomography. BMC Bioinforma. 16, 367 (2015).Article. Google Scholar.33. Wang, Z. et al. Image quality assessment: from error visibility to structural similarity. IEEE Trans. Image Process. 13, 600–612 (2004).ADS.Article. Google Scholar.34. Schäfer, J., Lee, S. C. & Kienle, A. Calculation of the near fields for the scattering of electromagnetic waves by multiple infinite cylinders at perpendicular incidence. J. Quant. Spectrosc. Radiat. Transf. 113, 2113–2123 (2012).ADS.Article. Google Scholar.35. D’Agostino, S. et al. Enhanced fluorescence by metal nanospheres on metal substrates. Opt. Lett. 34, 2381–2383 (2009).ADS.Article. Google Scholar.36. Yurkin, M. A. et al. Discrete dipole simulations of light scattering by blood cells. (Universiteit van Amsterdam, 2007).37. Kuchel, P. W. & Fackerell, E. D. Parametricequation representation of biconcave erythrocytes. Bull. Math. Biol. 61, 209–220 (1999).Article. Google Scholar.38. Bhattacharya, D. & Sharma, A. Split step nonparaxial finite difference method for 3D scalar wave propagation. Opt. Quantum Electron. 39, 865–876 (2007).Article. Google Scholar. Download references .Author information.Affiliations.Ecole Polytechnique Fédérale de Lausanne, Optics Laboratory, CH1015, Lausanne, Switzerland.Joowon Lim., Ahmed B. Ayoub., Elizabeth E. Antoine. & Demetri Psaltis.Authors.Search for Joowon Lim in:.PubMed • . Google Scholar .Search for Ahmed B. Ayoub in:.PubMed • . Google Scholar .Search for Elizabeth E. Antoine in:.PubMed • . Google Scholar .Search for Demetri Psaltis in:.PubMed • . Google Scholar .Contributions.J.L. carried out the algorithm modeling and computations. A.A. built the optical setup and carried out the optical experiments. E.A. prepared the samples. D.P. supervised the project. All authors contributed to the discussion and wrote the paper.Corresponding author.Correspondence to Joowon Lim.Ethics declarations. Conflict of interest. The authors declare that they have no conflict of interest.Supplementary information. Supplementary Information..Rights and permissions. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.Reprints and Permissions.About this article.Received.12 February 2019Revised.14 August 2019Accepted.19 August 2019Published.11 September 2019DOI.https://doi.org/10.1038/s4137701901951.

From：


