Assessment of learning tomography using Mie theory

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.


I. INTRODUCTION
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 [16,18,19] 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 [20].
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 [21]. 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.

A. Optical diffraction tomography (ODT)
The fundamental Helmholtz equation describing the scattering in inhomogeneous medium can be written as where U(r) is the total electric field: the sum of the incident field U i (r) and the scattered field U s (r). F (r) = k 2 4π n(r) 2 n 2 0 − 1 is the scattering potential in a sample immersed in a medium with refractive index n 0 . The optical wavelength in free space is λ resulting in wavenumber, k = 2πn 0 λ . The integral solution of Eq. (1) can be obtained using the homogeneous Green's function resulting in where G(r − r ′ ) = e ik|r−r ′ | |r−r ′ | is the Green's function of 3-D Helmholtz equation, Eq. 1.

B. Linear forward model
The scattered field U s (r) in Eq.
(2) is linear in U i (r) but nonlinear in F (r). Born and Rytov are the two linearizations. Both of them are based on the Wolf transform [11] as, where U Approx s (r) is determined depending on the approximation (Born : U Approx ). Under the assumption of planar incidence, U i (r) = e ik in ·r , Eq. 3 can be transformed to the following, where U Approx s (r; z = 0) is the measurement in the image plane and k z = k 2 − k 2 x − k 2 y . Eq. (4) 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.

C. 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(r) for U(x, y, z) = A(x, y, z)e ikz , BPM can be written as where ∆n is the contrast between RI of the sample (n(x, y, z)) and n 0 ,

D. 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 [18]. Specifically, the cost function is defined as where y ∈ C M denotes the experimental measurements, f ∈ R N denotes the object function, A : R N → C M is the forward model which can be either linear or nonlinear, (∇ x f ) 2 + (∇ y f ) 2 + (∇ z f ) 2 (∇ x , ∇ y , and ∇ z are finite difference operators in x, y and z dimensions, respectively.), N(f ) is the indicator function, N(f ) = ∞ (f < 0) or N(f ) = 0 (f ≥ 0), and τ 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) [19]. For LT, we also rely on the the stochastic gradient method to speed up the computation.

A. Simulation setup
To obtain the equivalent experimental measurements, we used the Mie theory to derive the scattered field by a single sphere [20,21]. The sample was illuminated at 95 different angles and the k-space trajectory of illumination angles is shown in Fig. 1 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, where f non and f lin are solutions acquired from the nonlinear and linear models, respectively.

C. 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 (∆x = ∆y = ∆z) of 0.08µm (single bead and multiple cylinders) or 0.0856µm (two beads). The regularization parameter, τ , was manually set as described in Table 1 (∆n : RI contrast, L : diameter in µm, ∆n 0 : 0.0406, and l 0 : 5 µm). For both models, the step size was reduced after every iteration, γ k+1 = 0.985γ k (single sphere and two beads) or γ k+1 = 0.99γ 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.

A. 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, Errors are in the appendix). The relative effect of the increase in RI contrast can be investigated by looking at Relative Error ( Fig. 3(b)). As the RI contrast increases, Relative Error increases in the range from 0.5 π to 1.5 π.

B. 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, Errors are in the supplementary materials). The advantage of LT against linear tomography can be clearly seen in Fig. 5  converging points of each Relative Error, it can be observed that the incremental change in Relative Error gets bigger with the increase of the diameter.

C. 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 sampleinduced phase delay to π. 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, Errors are in the appendix). In Case 3 where the two factors have compensatory effects on Errors, it was possible to observe that Relative Error values for the cases of 2.5 µm, 3.75 µm and 5 µm converge ( Fig. 7(b)) with the number of iteration.

D. Comparative analysis
The converging points of Errors and Relative Errors are plotted in Fig. 8. Fig. 8(a,c,e) not only show overall tendencies but also give us more information about how Errors converge depending on each individual factor explaining changes in Relative Errors as shown in Fig. 8(b,d,f). Overall, linear model breaks down with either the increase of the RI con- trast 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

E. Multiple cylinders
Even though each cylinder is a week scatterer, when they are aggregated, multiple scattering caused among them can introduce severe distortions [14]. To test the capacity for handling multiple scattering, simulated measurements for multiple cylinders has been generated using the Mie theory [21]. A set of three cylinders whose diameter is 4 µ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 Errors and Relative Errors, as shown in Fig. 10. Relative Errors for multiple cylinders rotated from 0 degree to 90 degrees.

F. 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 π 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 π. 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 [22]. 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).

G. Experimental results
We conducted an experiment to assess LT. Two 4.45 µ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 π, phase unwrapping must be deployed in the Rytov algorithm. When the phase unwrapping algorithm fails as the optical path across the 5 µm sphere increased to 3 π, 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.