2.1.

Experimental Setup

Figure 1 shows the schematic diagram of the experimental setup. In our system,²³ the sample placed between two cover glasses is illuminated sequentially at multiple angles and the scattered light is holographically recorded. A laser beam ($\lambda =561\mathrm{nm}$) is split into sample and reference arms by the first beam splitter. In the sample arm, a galvo mirror varies the angle of illumination on the sample using the 4F system created by L1 and OB1. The light transmitted through the sample is imaged onto the CMOS camera via the 4F system created by OB2 and L2. The beam splitter (BS2) recombines the sample and reference laser beams, forming a hologram at the image plane. The numerical apertures (NAs) of OB1 and OB2 are 1.45 and 1.4, respectively. For data acquisition, we capture multiple tomographic phase images by near-plane-wave illumination (Gaussian beam) with equally spaced incident angles. We use a differential measurement between the phase on a portion of the field of view on the detector that does not include the cell and the cell itself to maintain phase stability. Accordingly, complex amplitudes extracted from the measurements constitute the training set of our proposed CNN.

Fig. 1

Experimental setup (BS, beam splitter; GM, galva mirror; L, lens; M, mirror; and OB, objective) and measured hologram by CMOS. Scale bar, $10\mu \mathrm{m}$.

2.2.

Beam Propagation Method

We build the CNN, based on the forward model of light propagation,^21,23 to model the diffraction and propagation effects of light-waves. It is known that the scalar inhomogeneous Helmholtz equation completely characterizes the light field at all spatial positions in a time-independent form

2.3.

GPU-Based Implementation of CNN

A schematic architecture of our CNN is shown in Fig. 2. For constructing our neural network, we divide the computational sample space into thin slices with the sampling interval $\delta z$ along the propagation direction $z$. One slice corresponds to one layer in CNN. Within each layer, neurons specify the discretized light-field with transverse sampling intervals $\delta x$ and $\delta y$, respectively. The input layer is the incident field upon the sample. In terms of the Eq. (3), inputs are then passed from nodes of each layer to the next, with adjacent layers connected by alternating operations of convolution and multiplication. At the very last layer of our CNN, the output complex field amplitude is then bandlimited by the NA of the imaging system composed of lenses OB2 and L2 in Fig. 1. We implement the proposed network on the basis of TensorFlow framework. The connection weight $\delta n(\mathbf{r})$ can be trained using the Adam algorithm for optimization on the following minimization problem:

Fig. 2

Detailed schematic of our CNN architecture, indicating the number of layers ($Nz$), nodes ($Nx\times Ny$) in each layer and operations between adjacent layers. Here, $\ker ({\omega }_x,\hspace{0.17em}{\omega }_y)$ signifies $({\omega }_x^2+{\omega }_y^2)/(k_0{n}_0+\sqrt{k_0^2{n}_0^2-{\omega }_x^2-{\omega }_y^2})$ and we take the $\delta n(\mathbf{r})$ of a polystyrene bead as example.

3.1.

Tomographic Reconstruction of Simulated Data

In simulation, we consider a situation of three $5\mu \mathrm{m}$ beads of refractive index $n=1.548$ immersed into oil of refractive index $n_0=1.518$ shown in Fig. 3. The centers of the beads are placed at $(0,0,-3)$, (0, 0, 3), and (0, 5, 0), respectively, with the unit of micron. The training set of the framework is simulated as 81 complex amplitudes extracted from the digital-holography measurements with different angles of incidence evenly distributed in $[-\pi /4,\pi /4]$ by BPM, whereas the illumination is tilted perpendicular to the $x$-axis and the angle is specified with respect to the optical axis $z$. The size of the reconstructed volume is $23.04\mu \mathrm{m}\times 23.04\mu \mathrm{m}\times 23.04\mu \mathrm{m}$, with the sampling steps of $\delta x=\delta y=\delta z=144\mathrm{nm}$. For the network hyperparameters, we choose 600 training iterations in our GPU-based implementation with the batch size of 20, the initial learning rate of 0.001, and the regularization coefficient of $\tau =1.5$. The reconstructed results by our method and other reconstruction methods are shown in Fig. 4. The SNR defined in Ref. 21 of our result is 25.56 dB, 14 dB higher than the previous works. We can also observe much sharper edges of the reconstructed beads at the interface with less noise in the background from Fig. 5. The comparison between the proposed loss function and other regularized loss functions proves the higher reconstruction quality of the ${\ell }_1$-norm constraint than the ${\ell }_2$-case directly.

Fig. 3

Simulation geometry comprising three spherical beads with a refractive index difference of 0.03 compared with the background.

Fig. 4

Reconstruction results of three $5\mu \mathrm{m}$ beads. Comparison of the cross-sectional slices of the 3-D refractive index distribution of the sample along the $x-y$, $x-z$, and $y-z$ planes reconstructed by (a) proposed CNN, (b) CNN with ${\ell }_1$ fitting, ${\ell }_2$ regularization (L1 L2) and the regularization coefficient of 5, (c) CNN with ${\ell }_2$ fitting, ${\ell }_1$ regularization (L2 L1) and the regularization coefficient of 0.1, (d) learning approach²³ implemented on the same CNN settings (LA) with the regularization coefficient of 0.6, (e) optical diffraction tomography based on the Rytov approximation (ODT)¹³ with the positivity constraint and 100 iterations, and (f) iterative reconstruction based on the filtered backprojection method⁴ with the positivity constraint and 400 iterations. Scale bar, $5\mu \mathrm{m}$.

Fig. 5

Comparison of the refractive index profiles along the $z$-axis reconstructed by different algorithms and the ground truth.

In addition, we analyze the performance of our method under different noise levels and reduced sampling angles. For the noise test, we add Gaussian noise of different power levels to the 81 simulated measured complex amplitudes, which are represented as different SNRs of the training data. From the curve of the reconstructed SNR versus the noise level, as shown in Fig. 6(a), we can find our method maintains more robustness to the noise than other methods. This is especially useful in the case of shorter exposure time for higher scanning speed, where the data are always readout-noise limited. For the test of reduced sampling angles, we keep the range of incident angles fixed from $-\pi /4$ to $\pi /4$. The total number of the incident angles for the network training decreases from 81. The curve of the reconstructed SNR versus the number of the incident angles is shown in Fig. 6(b). Even with as few as 11 incident angles, we can still achieve comparable performance as the previous methods with 81 angles. This nearly eight-time improvement facilitates the development of high-speed 3-D refractive index imaging.

Fig. 6

Performance analysis for proposed approach with the same hyperparameter selection. (a) The curve of the reconstructed SNR versus the noise level and (b) the curve of the reconstructed SNR versus the number of the incident angles.

3.2.

Tomographic Reconstruction of a Biological Sample

To further validate the capability of the network, we display the experimental results on HeLa cells performed by our tomographic phase microscope as shown in Fig. 1. In detail, we illuminate the HeLa cells in culture medium of refractive index $n_0=1.33$ from 41 incident angles evenly distributed from $-35\mathrm{deg}$ to 35 deg. The measured hologram with an incident angle of 0 deg is shown in Fig. 1. The reconstructed volume is $36.86\mu \mathrm{m}\times 36.86\mu \mathrm{m}\times 23.04\mu \mathrm{m}$, composed of $256\times 256\times 160$ voxels (with the voxel size of $144\mathrm{nm}\times 144\mathrm{nm}\times 144\mathrm{nm}$). After the selection of hyperparameters, we set the regularization coefficient to $\tau =5$, with training iterations of 100, the batch size of 20, and the initial learning rate of 0.002. The performance comparison of different methods under different levels of data reduction is shown in Fig. 7. More details can be observed by our method even with fewer incident angles. Moreover, many fewer artifacts and noises exist in our results with large data reduction than other methods, which can be seen apparently in Fig. 8.

Fig. 7

Comparison of three reconstruction algorithms for various levels of data reduction on a HeLa cell. (a–d) Proposed CNN, (e–h) LA with the regularization coefficient of 1.5, (i–l) ODT with the positivity constraint and 20 iterations, (a, e, and i) 41 training data, (b, f, and j) 21 training data, (c, g, and k) 11 training data, and (d, h, and l) 6 training data. Scale bar, $10\mu \mathrm{m}$.

Fig. 8

Comparison of the reconstructed HeLa cell refractive index profiles along the $z$-axis with 41 training data.

4.1.

Comparison between $L_1$-Norm and $L_2$-Norm for Loss Function Design

Loss function design is a crucial element of learning algorithms, which determines the training process of the neural network. Regularized loss function comprises one data-fidelity term and one regularization term.

To the best of our knowledge, the presented study is the first to employ ${\ell }_1$ fitting for the regularized TPM. Generally, the choice of the data-fidelity term depends on the specified noise distribution. However, it is particularly common for solving normal image restoration problems, as under various constraints, images are always degraded with mixed noise and it is impossible to identify what type of noise is involved. The ${\ell }_2$-norm fitting relies on strong assumptions on the distribution of noise: there are no heavy tails and the distribution is symmetric. If either of these assumptions fails, then the use of ${\ell }_2$-norm is not an optimal choice. On the other hand, for the so-called robust formulation based on ${\ell }_1$-norm fitting, it has been shown that the corresponding statistics can tolerate up to 50% false observations and other inconsistencies.²⁷ Hence, ${\ell }_1$-norm data-fidelity term relaxes the underlying requirements for the ${\ell }_2$-case and is well suited to biomedical imaging especially when the noise model is undetermined [as shown in Figs. 4(a)–4(d)] and mixed [as shown in Fig. 6(a)].

As for the regularization term, we finally choose the anisotropic total variation (TV) regularizer in our method, which is an ${\ell }_1$ penalty directly on the image gradient. It is a very strong regularizer, which offers improvements on reconstruction quality to a great extent compared with the isotropic counterpart (${\ell }_2$ penalty).²⁹ Therefore, the edges are better preserved, which can be seen apparently from the comparison between Figs. 4(a) and 4(b).

4.2.

Selection of Hyperparameters

Selection of hyperparameters has an important effect on tomographic reconstruction results. In practice, many learning algorithms involve hyperparameters (10 or more), such as initial learning rate, minibatch size, and regularization coefficient. Reference 31 introduces a large number of recommendations for training feed-forward neural networks and choosing the multiple hyperparameters, which can make a substantial difference (in terms of speed, ease of implementation, and accuracy) when it comes to putting algorithms to work on real problems. Unfortunately, optimal selection of hyperparameters is challenging due to the high computational cost when using traditional regularized iterative algorithms.^21,23

In this study, our GPU-based implementation of CNN runs computation-intensive simulations at low cost and is possible to make the optimization of hyperparameters a more reproducible and automated process with modern computing facilities. Thus, we can gain better and more robust reconstruction performance with the GPU-based learning method. During the simulation and experiment, selection of hyperparameters varies with the biological sample and the range of incident angles. To make a convincing comparison, optimal hyperparameters have been selected for all the other reconstruction methods. The refractive index difference $\delta n(\mathbf{r})$ is initialized with a constant value of 0 for all the methods, and different optimal regularization coefficients are chosen for different regularized loss functions due to the different combinations of data-fidelity term and regularization term. The number of iterations is set to guarantee the convergence of each method, as shown in Fig. 9.

Fig. 9

Reconstructed SNR plotted as a function of the number of iterations for different reconstruction methods on simulated data. Hyperparameters are declared in Sec. 3.1.

4.3.

Subgradient Method and Adam Algorithm

In convex analysis,³⁰ the subgradient generalizes the derivative to functions that are not differentiable. A vector $g\in {\mathbb{R}}^n$ is a subgradient of a convex function at $x$ if

If $f$ is convex and differentiable, then its gradient at $x$ is a subgradient. But a subgradient can exist even when $f$ is not differentiable at $x$. There can be more than one subgradient of a function $f$ at a point $x$. The set of all subgradients at $x$ is called the subdifferential, and is denoted by $\partial f(x)$. Considering the absolute value function $|x|$, the subdifferential is

Subgradient methods are subgradient-based iterative methods for solving nondifferentiable convex minimization problems.

Adam is an algorithm for first-order (sub)gradient-based optimization of stochastic objective functions, based on adaptive estimates of lower-order moments. The method is aimed toward machine learning problems with large datasets and/or high-dimensional parameter spaces. The method is also appropriate for nonstationary objectives and problems with very noisy and/or sparse gradients. Adam works well in practice and compares favorably with other stochastic optimization methods regarding the computational performance and convergence rate.²⁶ It is straightforward to implement, is computationally efficient, and has little memory requirements, which is robust and well suited to TPM. Compared with the stochastic proximal gradient descent (SPGD) algorithm reported in Ref. 21, our GPU-based CNN trained with the Adam algorithm for optimization converges to the same SNR level and achieves twice the rate of convergence as shown in Fig. 10. To show the higher convergence rate of Adam fairly, we use the proposed ${\ell }_1$-norm loss function for both the Adam and SPGD training processes here, thus producing the same reconstructed SNR after convergence.

Fig. 10

Reconstructed SNR of proposed approach plotted as a function of the number of iterations for two different training optimization algorithms on simulated data with the same hyperparameters declared in Sec. 3.1.