Assessment of learning tomography using Mie theory

TL;DR

This paper evaluates a nonlinear reconstruction method called Learning Tomography (LT) in optical diffraction tomography, demonstrating its improved accuracy and robustness over linear models like Rytov, especially in complex multiple scattering scenarios.

Contribution

The paper introduces a simulation-based assessment of Learning Tomography using Mie theory, showing its advantages over traditional linear methods in handling multiple scattering and distortions.

Findings

LT outperforms linear models in accuracy and robustness.

LT corrects phase unwrapping distortions.

LT effectively handles multiple scattering in experiments.

Abstract

In Optical diffraction tomography, the multiply scattered field is a nonlinear function of the refractive index of the object. The Rytov method is a linear approximation of the forward model, and is commonly used to reconstruct images. Recently, we introduced a reconstruction method based on the Beam Propagation Method (BPM) that takes the nonlinearity into account. We refer to this method as Learning Tomography (LT). In this paper, we carry out simulations in order to assess the performance of LT over the linear iterative method. Each algorithm has been rigorously assessed for spherical objects, with synthetic data generated using the Mie theory. By varying the RI contrast and the size of the objects, we show that the LT reconstruction is more accurate and robust than the reconstruction based on the linear model. In addition, we show that LT is able to correct distortion that is evident in Rytov approximation due to limitations in phase unwrapping. More importantly, the capacity of LT in handling multiple scattering problem are demonstrated by simulations of multiple cylinders using the Mie theory and confirmed by experimental results of two spheres.

Tables

Table 1. Table 1: Parameters

Data (single type) Size (X$\times$Y$\times$Z)	Parameter	Linear	Nonlinear
Single sphere (350$\times$350$\times$128)	$\gamma$	0.01	0.001
	$\tau$	$\tau =2\Delta n\times L/\left(\Delta n_0\times l_0\right)$	0.1
Multiple cylinders (1024$\times$256)	$\gamma$	0.01	0.001
	$\tau$	0.4	0.1
Two beads (256$\times$256$\times$256)	$\gamma$	0.01	0.001
	$\tau$	4	1.5

Equations

$${\nabla }^2{U}_s(\text{r})+k^2{U}_s(\text{r})=-4\pi F(\text{r})U(\text{r}),$$

$$U_s(\text{r})={\int }_{V}F(\text{r′})U(\text{r′})G(\text{r}-\text{r′})d\text{r′},$$

$$U_s^{Approx}(\text{r})={\int }_{V}F(\text{r′})U_i(\text{r′})G(\text{r-r′})d\text{r′},$$

$$\frac{k_z}{2\pi i}{\iint }_{-\infty }^{+\infty }{U}_s^{Approx}(\text{r};z=0)e^{i(k_{x}x+k_{y}y)}dxdy$$

$$={\int }_{V}F(\text{r′})e^{-i(\text{k}-\text{kin})\cdot \text{r′}}d\text{r′},$$

$$A(x,y,z+dz)$$

$$\times F_{2D}^H[F_{2D}\{$$

$$C(f)=\frac{1}{2}||y-A(f)|{|}_2^2+\tau D(f)+N(f),$$

$$Error(f_{recon},f_{true})=\frac{||f_{recon}-f_{true}|{|}_2^2}{||f_{true}|{|}_2^2},$$

$$RelativeError(f_{lin},f_{non},f_{true})=\frac{Error(f_{non},f_{true})}{Error(f_{lin},f_{true})},$$

Full text

I Introduction

Quantitative phase imaging (QPI) is a tool for measuring the sample-induced phase delay, which relates to the refractive index (RI) contrast. Each material has its own distinct RI value and therefore QPI can provide physiological information Lee et al. (2013) such as the structure and dynamics of cells Popescu et al. (2006); Park et al. (2011), quantification of specific molecules Park et al. (2009); Jung et al. (2014) and the dry mass Barer (1952); Popescu et al. (2008), which serves as an important factor for monitoring cell growth. Among other QPI techniques, optical diffraction tomography (ODT) enables us to visualize 3-D RI distributions from multiple 2-D scattered fields acquired at various illumination angles. It provides the physiological information through measuring the 3-D RI distribution without any exogenous labeling agents, making it a powerful tool for various physiological studies.

Despite the many potential applications of ODT, there are obstacles that need to be resolved. One obstacle comes from the fact that views far from the optical axis are usually not accessible due to the limited numerical aperture of the optics. This is referred to as the missing cone problem. It causes the underestimation of RI values and the elongation of RI tomograms along the optical axis Sung and Dasari (2011); Sung et al. (2009). The missing cone problem can be mediated using sparsity-based regularization algorithms Sung and Dasari (2011); Lim et al. (2015). A more fundamental problem is the inherent nonlinearity in scattering through inhomogeneous media Born and Wolf (1999). These nonlinear effects are assumed to be negligible in conventional linear ODT frameworks (Born and Rytov approximations). The validity of the approximations were well studied in showing the restriction of the linear assumption Slaney et al. (1984); Chen and Stamnes (1998), in simple objects such as a single cell. In addition, distortions caused by multiple scattering can be severe for thick or high contrast samples especially when several scattering objects are aggregated Azimi and Kak (1983).

Recently, reconstruction algorithms that consider the nonlinear process have been proposed Tian and Waller (2015); Kamilov et al. (2015). The Beam propagation method (BPM) can be used as the forward model combined with the sparsity based regularization in the iterative reconstruction scheme Kamilov et al. (2015). BPM is a method that describes the propagation of light through a medium. It consists of two sub-steps: Diffraction followed by a refraction step Okamoto (2010). The BPM can implement multiple scattering at different depths within the medium. Therefore, iterative reconstruction algorithms that combine the BPM and a sparsity-based regularization can outperform the conventional linear model (Rytov approximation) Kamilov et al. (2015, 2016).

However, to fairly investigate the capacity of the nonlinear model over that of the linear model, the same reconstruction framework should be used except for the forward model part. In order to distinguish the performance of LT from the effects of a sparsity-based regularization, it should be compared with the linear tomography which utilizes the same reconstruction algorithm scheme including the sparsity-based regularization. In addition, in order to quantitatively compare any improvements made by the nonlinear algorithm, we need the ground truth of the 3-D RI distribution of the object.

In this paper, the same algorithmic scheme Kamilov et al. (2015, 2016); Beck and Teboulle (2009) was used for the linear and nonlinear algorithms except for the forward model part. Simulated measurements were generated using the Mie theory whose analytical solution serves as the ground truth Bohren and Huffman (2008). The two models are tested on various schemes where the RI contrast and the size of the sample are controlled under three different schemes: 1) variable RI contrast with fixed size, 2) variable size with fixed RI contrast, and 3) different RI contrast and size with fixed sample-induced phase delay. In addition, in order to compare the ability of each model in dealing with multiple scattering caused by multiple objects, the two models are tested using the Mie theory for multiple cylinders Schäfer et al. (2012). After that, the capacity of each model when phase unwrapping fails is discussed. Finally, we applied algorithms on an experimental data confirming the simulation results.

II Theory

II.1 Optical diffraction tomography (ODT)

The fundamental Helmholtz equation describing the scattering in inhomogeneous medium can be written as

[TABLE]

where $U(\textbf{r})$ is the total electric field: the sum of the incident field $U_{i}(\textbf{r})$ and the scattered field $U_{s}(\textbf{r})$. $F(\textbf{r})=\frac{k^{2}}{4\pi}\left(\frac{n(\textbf{r})^{2}}{n_{0}^{2}}-1\right)$ is the scattering potential in a sample immersed in a medium with refractive index $n_{0}$. The optical wavelength in free space is $\lambda$ resulting in wavenumber, $k=\frac{2\pi n_{0}}{\lambda}$ . The integral solution of Eq. (1) can be obtained using the homogeneous Green’s function resulting in

[TABLE]

where G(\textbf{r}-\textbf{r{}^{\prime}})=\frac{e^{ik|\textbf{r}-\textbf{r{}^{\prime}}|}}{|\textbf{r}-\textbf{r{}^{\prime}}|} is the Green’s function of 3-D Helmholtz equation, Eq. 1.

II.2 Linear forward model

The scattered field $U_{s}(\textbf{r})$ in Eq. (2) is linear in $U_{i}(\textbf{r})$ but nonlinear in $F(\textbf{r})$. Born and Rytov are the two linearizations. Both of them are based on the Wolf transform Born and Wolf (1999) as,

[TABLE]

where $U_{s}^{Approx}(\textbf{r})$ is determined depending on the approximation (Born : $U_{s}^{Approx}(\textbf{r})$ = $U(\textbf{r})-U_{i}(\textbf{r})$ , Rytov : $U_{s}^{Approx}(\textbf{r})$ = $U_{i}(\textbf{r})\log\frac{U(\textbf{r})}{U_{i}(\textbf{r})}$). Under the assumption of planar incidence, U_{i}(\textbf{r})=e^{i\textbf{k{}^{in}}\cdot\textbf{r}}, Eq. 3 can be transformed to the following,

[TABLE]

where $U_{s}^{Approx}(\textbf{r};z=0)$ is the measurement in the image plane and $k_{z}=\sqrt{k^{2}-k_{x}^{2}-k_{y}^{2}}$. Eq. (II.2) is equivalent to the 3-D k-space of the object filled with the 2-D Fourier transform of the measured field over the Ewald’s sphere.

II.3 Nonlinear forward model

We mow describe the BPM based reconstruction algorithm for optical diffraction tomography. We can propagate the light through inhomogeneous medium by splitting the process in multiple fine steps, where each step consists of a diffraction step followed by a refraction step. Denoting the envelope of the wave as $A(\textbf{r})$ for $U(x,y,z)=A(x,y,z)e^{ikz}$, BPM can be written as

[TABLE]

where $\Delta n$ is the contrast between RI of the sample ($n(x,y,z)$) and $n_{0}$, $k_{z}=\sqrt{k^{2}-k_{x}^{2}-k_{y}^{2}}$. The 2-D Fourier transform operator ($F_{2D}$) is defined as $A(k_{x},k_{y},z)=\iint A(x,y,z)e^{-i(k_{x}x+k_{y}y)}dxdy$ .

II.4 Iterative reconstruction algorithms

Once the forward model is determined as either linear or nonlinear, we can specify a cost function that combines an error term and a regularization that incorporates prior knowledge about the sample. We impose the total variation (TV) and the non-negativity constraints Kamilov et al. (2016). Specifically, the cost function is defined as

[TABLE]

where $y\in\mathbb{C}^{M}$ denotes the experimental measurements, $f\in\mathbb{R}^{N}$ denotes the object function, $A:\mathbb{R}^{N}\rightarrow\mathbb{C}^{M}$ is the forward model which can be either linear or nonlinear, $D(f)=\displaystyle\sum_{n=1}^{N}\sqrt{(\nabla_{x}f)^{2}+(\nabla_{y}f)^{2}+(\nabla_{z}f)^{2}}$ ($\nabla_{x},\nabla_{y},\text{ and }\nabla_{z}$ are finite difference operators in $x,y\text{ and }z$ dimensions, respectively.), $N(f)$ is the indicator function, $N(f)=\infty$ $(f<0)$ or $N(f)=0$ $(f\geq 0)$, and $\tau$ is the regularization parameter, setting the relative weight of the regularization term. To minimize the cost function, Eq. 6. In both cases, we used the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) Beck and Teboulle (2009). For LT, we also rely on the the stochastic gradient method to speed up the computation.

III Method

III.1 Simulation setup

To obtain the equivalent experimental measurements, we used the Mie theory to derive the scattered field by a single sphere Bohren and Huffman (2008); Schäfer et al. (2012). The sample was illuminated at 95 different angles and the k-space trajectory of illumination angles is shown in Fig. 1. To investigate the performance of each forward model, we tested the algorithms on following three different schemes. The first is the case of different RI contrasts with fixed radius (Case 1), the second is the case of different radius with fixed RI contrast (Case 2). The third is the case with a fixed sample-induced phase delay but different in both RI contrast and diameter (Case 3). In addition, to investigate the capacity of LT in handling multiple scattering caused by multiple samples, Mie theory for multiple cylinders was used Schäfer et al. (2012) to investigate both forward models. Uniformly distributed 101 angles between $-\pi/4$ and $\pi/4$ were used to scan the sample. To quantitatively evaluate the performance, we calculated the error defined as,

[TABLE]

where $f_{recon}$ is the reconstructed solution and $f_{true}$ is the ground truth from the Mie simulation. In addition, to compare the relative performances of two forward models, the relative error is defined as,

[TABLE]

where $f_{non}$ and $f_{lin}$ are solutions acquired from the nonlinear and linear models, respectively.

III.2 Experimental setup

The experimental apparatus consists of a Mach-Zehnder interferometer. The signal and reference arm are recombined at an angle of 1.43 degrees before the detector in order to record off-axis digital holograms. The light source is a continuous wave laser diode at 406 $nm$ with a coherence length of 250 $\mu m$. The light was spatially filtered and collimated. The sample consisted of two polystyrene microspheres with a nominal diameter of 4.45 $\mu m$. The spheres were placed on two 150 $\mu m$ thick glass cover slips facing each other, each sphere sitting on a different cover slip. The gap between the cover slips was filled with a low fluorescence immersion oil (Nikon type N, n0 = 1.518 at 546.1 $nm$ and ve = 41), with a refractive index of 1.5370 at 406 $nm$. The coverslips and spheres assembly was placed in the signal arm between two infinity-corrected 100x oil-immersion microscope objectives (UPlanApo NA1.4 on the detection side and UplanFI NA1.3 on the illumination side). The same oil was used between the objectives and the cover slips. The effective numerical aperture of the system is 1.3. The image of the spheres was projected on a S-CMOS camera (Andor Neo) with an effective magnification of 111. The sphere were illuminated with plane wave through the illumination objective at different incidence angles. This could be achieved by placing two galvo-mirrors (one for the X and the other one for the Y direction) in conjugate image planes of the sample. We acquired 160 views equally spaced on a circle in the $k_{x}-k_{y}$ plane at an angle of 42 degrees from the optical axis (effective angle on the sample).

III.3 Reconstruction setup

The algorithms were implemented using custom scripts in MATLAB R2016b (MathWorks Inc., Natick, MA, USA) on a desktop computer (Intel Core i7-6700 CPU, 3.4 GHz, 32 GB RAM). For faster computation, a graphic processing unit (GPU, GeForce GTX 1070) was utilized. The computational space was sampled with the step size ($\Delta x=\Delta y=\Delta z$) of $0.08\mu m$ (single bead and multiple cylinders) or $0.0856\mu m$ (two beads). The regularization parameter, $\tau$, was manually set as described in Table 1 ($\Delta n$ : RI contrast, $L$ : diameter in $\mu m$, $\Delta n_{0}$ : 0.0406, and $l_{0}$ : 5 $\mu m$). For both models, the step size was reduced after every iteration, $\gamma^{k+1}=0.985\gamma^{k}$ (single sphere and two beads) or $\gamma^{k+1}=0.99\gamma^{k}$ (multiple cylinders), and the iteration numbers for TV and FISTA are 20 and 200, respectively. In case of LT, the stochastic gradient method was used with 8 randomly chosen angles out of total angles selected at each iteration. For every case, the first order Rytov approximation was used as the initial condition.

IV Results

IV.1 Case 1

Case 1 shows the effect of the RI contrast on the performance of the linear and nonlinear models. While both models show gradual distortions with the increase of the RI contrast in each case, the nonlinear model outperforms the linear one (Fig. 2). This is quantified by the $Error$ curve in Fig. 3(a), showing lower $Error$ of the nonlinear model for all cases (For the other cases, $Error$s are in the appendix). The relative effect of the increase in RI contrast can be investigated by looking at $Relative\text{ }Error$ (Fig. 3(b)). As the RI contrast increases, $Relative\text{ }Error$ increases in the range from 0.5 $\pi$ to 1.5 $\pi$.

IV.2 Case 2

For Case 2, we assessed the validity of the models by increasing the diameter of sphere. It can be clearly observed that the reconstruction results from the linear model tend to shrink along the optical axis with the increase of the diameter (Fig. 4). In contrast, the nonlinear model shows consistent tomograms (Fig. 4) and lower $Error$ values than that of the linear model (Fig. 5(a), for the other cases, $Error$s are in the supplementary materials). The advantage of LT against linear tomography can be clearly seen in Fig. 5(b). We can observe the gradual increase of $Relative\text{ }Error$ values with the gradual increase of the diameter, indicating the gradual break-down of the linear model. To be specific, by looking at the converging points of each $Relative\text{ }Error$, it can be observed that the incremental change in $Relative\text{ }Error$ gets bigger with the increase of the diameter.

IV.3 Case 3

The purpose of Case 3 is to see the relative importance of two different factors, the RI contrast and the size of the sample. Here, the size of sphere varies keeping the sample-induced phase delay to $\pi$. Therefore, each sample differs in both RI contrast and size, but the product of the two factors remains the same. It is obvious that LT outperforms linear tomography as shown in Fig. 6 and Fig. 7(a) (For the other cases, $Error$s are in the appendix). In Case 3 where the two factors have compensatory effects on $Error$s, it was possible to observe that $Relative\text{ }Error$ values for the cases of 2.5 $\mu m$, 3.75 $\mu m$ and 5 $\mu m$ converge (Fig. 7(b)) with the number of iteration.

IV.4 Comparative analysis

The converging points of $Error$s and $Relative\text{ }Error$s are plotted in Fig. 8. Fig. 8(a,c,e) not only show overall tendencies but also give us more information about how $Error$s converge depending on each individual factor explaining changes in $Relative\text{ }Error$s as shown in Fig. 8(b,d,f). Overall, linear model breaks down with either the increase of the RI contrast or the size, but the dependence on each individual factor shows different patterns as you can see from the blue curves in Fig. 8(a,c). Compared to that, LT shows much less changes in converging points (red curves in Fig. 8(a,c)). While it shows positive slope in Fig. 8(a), interestingly, the slope remains negative even when the size increases. These effects result in different dependencies of $Relative\text{ }Error$s as in 8(b,d). When both of the RI contrast and the size change (Fig. 8(e,f), it become more complex. Converging points of $Relative\text{ }Error$s almost stay constant in the range from 0.5 $l_{0}$ to 1 $l_{0}$, and increase fast after 1 $l_{0}$ as shown in Fig. 8(f). It comes from the different effectiveness of each factor as shown in Fig. 8(e). To be specific, in terms of $Relative\text{ }Error$s, when the RI contrast increases and the size decreases, the effect of each factor compensates the opposite effect of the other factor relatively showing small changes in the converging points of $Relative\text{ }Error$. However, the effect of the increase in size is stronger than the effect of the decrease in the RI contrast with the increase in the size over the size, $l_{0}$.

IV.5 Multiple cylinders

Even though each cylinder is a week scatterer, when they are aggregated, multiple scattering caused among them can introduce severe distortions Azimi and Kak (1983). To test the capacity for handling multiple scattering, simulated measurements for multiple cylinders has been generated using the Mie theory Schäfer et al. (2012). A set of three cylinders whose diameter is 4 $\mu m$ and RI value is $0.0508$ were considered. Then, the set of cylinders was rotated from 0 degree to 90 degrees resulting in the gradual increase of the multiple scattering effect. As shown in Fig. 9, we can clearly see that LT outperforms the linear model for all cases. While the linear tomography shows leakage between cylinders, the LT maintains a clear distinction of three cylinders. It can be quantitatively confirmed through the converging points of $Error$s and $Relative\text{ }Error$s, as shown in Fig. 10.

IV.6 Phase unwrapping

The Rytov algorithm makes use of calculated on the measured phase of each projection. If the optical thickness of the sample exceeds 2 $\pi$ then the measured phase must be unwrapped. Here, we compare the performance of these models when the phase unwrapping algorithm fails. The RI contrast is increased so that sample-induced phase delay resulted in 3 $\pi$. It is then difficult to properly unwrap the phase, even with a state-of-the-art algorithm such as a PUMA which was used throughout this paper Bioucas-Dias and Valadao (2007). Since the Rytov approximation requires unwrapped phases for the reconstruction, the failure in unwrapping directly relates to severe distortion in the reconstruction result of Rytov approximation, as shown in Fig. 11(a). It is critical not only to linear tomography, which uses the linear model based on Rytov approximation, but also to the LT because it is more likely for LT to fall in a local minimum. Compared to the linear tomography which only shows denoising effects (Fig. 11(b)), the LT is able to successfully reconstruct the tomograms (Fig. 11(c)) even though it uses 11(a) as the initialization. The effect is more dramatic in terms of $Error$ as in Fig. 11(d).

IV.7 Experimental results

We conducted an experiment to assess LT. Two 4.45 $\mu m$ were placed in a row so they overlapped in z-axis. We applied Linear tomography and LT on an experimental data. In case of linear tomography which is based on Rytov, a location of the image plane is very important. Fig. 9(c) confirms this fact. As the sample is placed farther from the image plane, RI tomograms gets either smaller in size and higher in RI contrast or bigger in size and lower in RI contrast. It was experimentally confirmed through 2 beads. It comes from the fact that Rytov directly uses phase perturbation in the field. As shown in Fig. 12(a-c), one below the image plane gets smaller in size and higher in RI contrast and the other above the image plane gets bigger in size and lower in RI contrast. By contrast, LT clearly reconstructs two beads of equal size and contrast (Fig. 12(d-f)). Since BPM describes the propagated field itself through the sample, it is not only unaffected by the location of the image plane but also able to handle multiple scattering.

V Conclusion

In this paper, we rigorously compare the LT against the conventional linear tomography. Mie theory provided the analytical solution for the scattered filed given a sphere so that we were able to evaluate the reconstruction fidelity of each model accurately. To investigate the capability of each model in dealing with nonlinearity, two factors which are directly related to nonlinearity, the RI contrast and the size, were controlled either independently or simultaneously. In all the cases, LT consistently outperformed linear reconstruction methods. We attribute this to the fact that the BPM used by LT captures multiple forward scattering events. The improvement in performance becomes more pronounced as the index or the diameter of the beads increases.  Our results indicate that LT provides a more pronounced improvement in the image quality for large size beads. The most dramatic improvement was observed when we imaged multiple objects (three cylinders). This confirms the observation that object size matters since we can consider  the set of 3 cylinders set as a single large object. In general the BPM performs relatively  well in the simulation of inhomogeneous media with small index contrast since the main limitation of the method is the assumption that reflections can be neglected. For samples whose optical path exceeds 2 $\pi$, phase unwrapping must be deployed in the Rytov algorithm. When the phase unwrapping algorithm fails as the optical path across the 5 $\mu m$ sphere increased to 3 $\pi$, the linear reconstruction becomes severely distorted. For this case, we observed that the iterations of the LT algorithm were able to correct the distortions that are evident in the Rytov reconstruction due to the phase unwrapping limitations. Finally, the reconstruction of the Rytov algorithm is in focus only at the plane of best focus of the optical system [?,our future paper]. For thick samples the sample becomes blurred away from the focal plane. This is evident in Figures 9(a-e) and 12(a-c). The distortion evident in Fig. 9(c) for example is a combination of two factors: Defocussing and multiple scattering. The BPM helps alleviate both of these problems since it takes multiple scattering into account  and it allows us to keep the entire sample in focus. It was experimentally confirmed using 2 Beads overlapping in the optical axis. As the sample becomes more complex (thicker, higher index contrast) ultimately the BPM provides an inadequate estimate for the scattered field by the object since reflections were neglected and the vectorial nature of the optical field was ignored. In this case the only way to realize an improvement in performance is by adopting a more sophisticated scattering model. There is an intermediate level of sample complexity however where BPM still provides a reasonably accurate prediction of the scattered field but the nonlinear inversion problem becomes very difficult due to emergence of strong local minima. In this regime we believe that there is a global minimum which is a good approximation of the true object but the LT algorithm cannot find it. It is possible that more powerful optimization algorithms than the stochastic decent we used in this paper can provide significant improvement in performance.

Appendix A Supplementary figures

As supplementary figures, $Error$s for each sample for each case are provided in detail.