1. Introduction

As a label-free, non-destructive and quantitative imaging technique, optical diffraction tomography (ODT) has been widely used to reconstruct the refractive-index (RI) distribution of phase objects, such as biological cells and tissues, and optical fibers [1–5]. ODT reconstructs the RI distributions of transparent samples by interferometrically measuring diffracted fields over multiple illumination angles and solving the inverse scattering problem, typically by adopting the weakly-scattering approximation [6–8]. When such an assumption does not hold, e.g., when the object has a high contrast, complicated structure or large optical path difference (OPD), conventional ODT inversion becomes ineffective [9,10]. To extend the application of ODT reconstruction to the multiply-scattering regime, optimization algorithms based on the beam propagation method [11–13] or the Lippmann-Schwinger equation were developed [14–17]; however, such techniques can be time-consuming or trapped in a local minimum, especially for objects with large OPDs. Recently, a computationally efficient iterative optical diffraction tomography (iODT) algorithm was proposed for the reconstruction of RI distributions of phase objects beyond the weakly-scattering regime [10].

While the reconstruction of RI distributions of transparent samples is useful for biological imaging, the absorption properties of tissue provide insight about light penetration and biological functions such as delineation of oxygenated and deoxygenated blood [18]. It was demonstrated that ODT can reconstruct objects with complex RI under the weakly-scattering assumption [19–23]. Since multiple scattering is common in biological samples, e.g., a cluster of cells can produce multiple scattering even if each cell is a weak scatterer [24], simultaneous reconstruction of the distributions of the refractive index and the absorption coefficient of multiply-scattering objects is highly desirable, albeit a hitherto unexplored topic, to the best of our knowledge.

In this paper, we will first explore the efficacy of iODT for reconstruction of objects with both multiple scattering and loss, i.e., complex RI. We have observed that for samples of higher RI contrast, or complicated structure, iODT reconstructions are susceptible to complex-valued artifacts, which may introduce “crosstalk” between the real and imaginary portions of reconstructions. To help mitigate the crosstalk error, we propose an error subtraction (ES) method to serve as an “add-on” module to iODT that extends the effectiveness of iODT, without changing the main skeleton of the algorithm. This new technique is hereafter referred to as iODT-ES.

2. Methods

The recently-proposed iODT algorithm (without the ES module) is described in Algorithm 1 and Fig. 1. The object function (viz. the scattering potential) is proportional to the permittivity difference between the object to be imaged and the background medium, i.e., $ f = {({2\pi {k_0}} )^2}({{n^2} - n_\textrm{b}^2} ),$ where ${k_0} = 1/{\lambda _0}$ is the wavenumber, ${\lambda _0}$ is the wavelength of the probe wave in free space, and ${n_\textrm{b}}$ is the background RI. iODT reconstructs the object function iteratively until the errors between the measured and computed fields are minimized. The true diffracted fields are measured experimentally [25] or obtained through simulation of phantoms using techniques such as the finite-difference time-domain (FDTD) method [26]. The global and local coordinates are related by $x = \xi \textrm{cos}\theta - \eta \textrm{sin}\theta ,{\; }y = \xi \textrm{sin}\theta + \eta \textrm{cos}\theta ,$ where $\theta $ is the illumination angle. The same symbol is used to represent both the function and its Fourier transform, with the argument specifying the domain. With the initial guess of the object function set to ${f^{(0 )}} = 0$, iODT reconstruction in the first iteration is identical to conventional ODT. Subsequently, the known incident field ${u_0}$ and the true diffracted fields ${u_\textrm{t}}$, shown in Fig. 1(a), are respectively forward and backward propagated through an estimate of the object using an efficient and accurate propagation solver ${\bf {\cal S}}$. As demonstrated in Fig. 1(b), the difference between these propagated fields within the entire reconstruction area produces the perturbative complex phase that forms a correction yielding the next estimate of the object function. The constraint ${\bf {\cal C}},$ e.g., non-negativity, can be applied to the reconstructed object if prior knowledge is available. Normalized root-mean-square (nRMS) errors between the simulated and true sinograms are defined as follows,

 

Fig. 1. Illustration of iODT. (a) The propagation module and the tomographic imaging configuration. $({x,\; y} )$: global coordinates; $({\xi ,\eta } )$: local coordinates; $\theta $: illumination angle; ${u_0}$: known incident wave; ${u_\textrm{t}}$: true diffracted fields measured experimentally or through simulations of phantoms; ${f^{(m )}}$: reconstructed object at the ${m^{\textrm{th}}}$ iteration; ${\bf {\cal S}}$: propagation solver. (b) Rytov-based inversion module. The forward and backward propagated fields produce perturbative complex phases, which are filtered in the frequency domain and summed in the spatial domain, to produce a perturbative correction to the object function.

Download Full Size | PDF

The true diffracted fields are either measured experimentally or computed through simulation of phantoms using FDTD method, which is well-known for its accuracy. As shown in Algorithm 1, or Fig. 1(a), the modular nature of iODT allows the user to define the propagation solver ${\bf {\cal S}}$, – as computational accuracy and efficiency considerations deem appropriate. To model the effects of multiple-scattering, we elect to use the wide-angle beam propagation method (WABPM) [27] for its computational efficiency and reasonable accuracy; however, it is less accurate than FDTD.

The discrepancy between the propagation solver ${\bf {\cal S}}$ and the solver producing the true diffracted fields, especially for samples of high-contrast or complicated structure, can result in errors in the forward-backward fields, and by extension the perturbative Rytov phase ${\Delta }{\varphi _\textrm{R}}$. According to the Rytov-based filtered backpropagation algorithm, as described in Algorithm 1, errors in the real and imaginary parts of ${\Delta }{\varphi _\textrm{R}}$ result in complex-valued errors in the perturbative object function correction, ${\Delta }f$. Consequently, the solver discrepancy (e.g., between WABPM and FDTD) can cause complex-valued artifacts in reconstruction, which in turn may introduce crosstalk artifacts between the reconstructed real and imaginary parts.

One way to suppress such crosstalk artifacts is to reduce the solver discrepancy, for example, by selecting a more accurate forward-backward solver, which more accurately models the true fields, computed from FDTD or measured in experiment. To demonstrate this, we create a complex-valued phantom whose real and imaginary parts are spatially separated, as shown in Fig. 2(a). Figure 2(b), the iODT reconstruction, where WABPM is used as the propagation solver ${\bf {\cal S}}$, is noticeably contaminated by real-to-imaginary crosstalk. We then repeat the experiment by changing the solver ${\bf {\cal S}}$ to FDTD, as plotted in Fig. 2(c), and observe a significant reduction in crosstalk contamination. Although the iODT reconstruction, using FDTD, is more accurate, FDTD is very time-consuming (about 200 times slower than WABPM), which suggests that such an approach for reduction solver error is less practical – especially for larger samples.

 

Fig. 2. Distributions of the real (upper row) and imaginary (lower row) parts of the object function of (a) phantom, (b) reconstruction using iODT with WABPM, (c) reconstruction using iODT with FDTD.

Download Full Size | PDF

An alternative approach for reduction of the crosstalk artifacts is to calibrate the error in iODT and then subtract it out. Since many samples, like soft tissue and optical fiber, are mostly transparent in the optical regime, the imaginary parts of the object functions of such samples are usually negligible, compared to the real parts. For this reason, the imaginary-to-real crosstalk is considered negligible compared to the real part of the object function itself. Exploring the reconstructed real part of the object function that is less sensitive to the crosstalk artifacts, we propose an ES method, serving as an add-on module to iODT, to help suppress the real-to-imaginary crosstalk artifacts caused by the solver discrepancy. This way, accurate reconstructions can be obtained without scarifying computational efficiency. The main idea of ES is that once a good estimate of the real part of the object function is obtained, as an approximation, the real-to-imaginary crosstalk induced by the real part can be subtracted out, as described in Algorithm 2 and Fig.  3. 

When iODT reconstruction (without the ES module) converges at the ${p^{\textrm{th}}}$ iteration, we assume that the real part of the reconstructed object $\textrm{Re}[{{f^{(p )}}} ]$ approximates the real part of the true object $\textrm{Re}[{{f_\textrm{t}}} ].$ However, the imaginary part $\textrm{Im}[{{f^{(p )}}} ]$ contains loss/gain artifacts. These artifacts occur because even if the object function is purely real, both ODT and iODT reconstructions produce an erroneous imaginary component in the object function, which is manifested as an artifact in the loss or gain of the object. In the ES module, the estimated phase object $\textrm{Re}[{{f^{(p )}}} ]$ without any loss/gain is plugged into iODT to compute the corresponding loss/gain artifact $g,$ which is then subtracted out for all subsequent iODT-ES iterations starting from ${f^{({p + 1} )}}$. To do so, fields $u_\textrm{t}^{\prime}$ diffracted by $\textrm{Re}[{{f^{(p )}}} ]$ are first calculated using FDTD for all illumination angles, and then iODT is employed to estimate the corresponding artifact $g$. In all subsequent iODT iterations ($m\;>\;p$), artifact g is subtracted from the imaginary component of the current estimate of the object function, and the real part of the reconstruction is reset to $\textrm{Re}[{{f^{(p )}}} ]$ in each iODT-ES iteration until the stopping condition is satisfied again.

3. Results

To test the effectiveness of the proposed iODT-ES method, we create multiply-scattering phantoms with complex RI distributions. The fields ${u_\textrm{t}}$, diffracted by the phantom over illumination angles from ${0^\textrm{o}}$ to ${355^\textrm{o}}$ with an increment of ${5^\textrm{o}}$, are calculated by the FDTD method. If the phantom to be imaged or the reconstructed object in each iteration has a rotational symmetry of order $\alpha $, only the diffracted fields corresponding to the angles from ${0^\textrm{o}}$ to ${360^\textrm{o}}/\alpha $ need to be computed. WABPM is employed as the propagation solver in iterative reconstruction. For conventional ODT reconstruction, we adopt the Rytov approximation, which has been shown to be superior to the Born approximation [9]. L2-norm phase unwrapping [28] is used to remove the phase ambiguities in the imaginary part of the perturbative complex phase $\Delta \varphi _\textrm{R}^{(m )}$ shown in Algorithm 1. The numerical validations of the reconstruction are here reported in terms of the complex RI ${n}$, because it is more intuitive, although the algorithm is formulated in terms of the equivalent object function ${f}$ (see Algorithm 2 or Fig. 3.) We choose the parameters $Q = 10$ and the error threshold ${\epsilon _{\textrm{tol}}} = {10^{ - 3}}$, respectively. All simulations are implemented in MATLAB 2018a using a computer equipped with an Intel Core i7-5820 K 3.3 GHz CPU and a 32GB RAM.

 

Fig. 3. Illustration of the iODT-ES algorithm. The true diffracted fields ${u_\textrm{t}}$ are measured experimentally or obtained through simulation of phantoms. With the initial guess ${f^{(p )}}$, the known incident field ${u_0}$ and the true fields ${u_\textrm{t}}$, iODT converges at the ${p^{\textrm{th}}}$ iteration and the reconstruction ${f^{(p )}}$ has a good estimate of the real part, while the real-to-imaginary crosstalk artifacts appear in the imaginary part. Based on the phase object $\textrm{Re}[{{f^{(p )}}} ]$ and its associated diffracted fields $u_\textrm{t}^{\prime}$ calculated by FDTD, iODT produces the reconstruction g, in which the crosstalk artifacts appear. With the phase object $\textrm{Re}[{{f^{(p )}}} ]$ as the estimate of the ${({p + 1} )^{\textrm{th}}}$ iteration and the true diffracted fields ${u_\textrm{t}}$, iODT produces the reconstruction ${f^{({p + 1} )}}$, where the crosstalk artifacts also occur. In the ES process, the real part of reconstruction is reset to be the phase object $\textrm{Re}[{{f^{(p )}}} ]$ obtained from iODT, while the imaginary part is subtracted by the estimated error $\textrm{Im}[g ]$. This ES process is repeated until the stopping condition is satisfied again.

Download Full Size | PDF

3.1 Shepp-Logan phantom

In the first example, a Shepp-Logan phantom with complex RI distribution is created as shown in Fig. 4(a). Although the real and imaginary parts of the RI distributions of physical objects usually have similar boundaries or orientations, we purposely make them different so that we can evaluate the reconstruction fidelity individually. The wavelength of the probe light is $\lambda_0 = 406$ nm. The RI of the background medium is ${n_\textrm{b}} = 1.333$, so that the wavelength in the background medium is $\lambda_\textrm{b}{ = \lambda_0}/{n_\textrm{b}} = 304.6$ nm. The real part of the RI has a contrast $({1.407 - {n_\textrm{b}}} )/{n_\textrm{b}} \approx 5.6\%$, and the imaginary part introduces approximately an average of 13.5% energy loss of the diffracted fields over all illumination angles. The reconstruction dimension is ${20\lambda_\textrm{b}} \times {20\lambda_\textrm{b}}$ with a grid size ${\Delta}x = {\Delta}y{ = \lambda_{\textrm{b}}}/20$. As shown in Fig. 4(b), conventional ODT reconstruction produces both real and imaginary parts with values smaller than the phantom as a result of the weakly-scattering assumption. Also, the crosstalk artifacts are apparent in the imaginary part. By contrast, the iODT result ($m = 28$) shown in Fig. 4(c) is superior to the ODT result in terms of the reconstructed values, shape, and artifact suppression. The data processing time of iODT reconstruction is about two minutes per iteration.

 

Fig. 4. Distributions of the real (upper row) and imaginary (lower row) parts of the complex RI of (a) Shepp-Logan phantom with an index contrast of 5.6%, (b) conventional ODT reconstruction, (c) iODT reconstruction at iteration $m = 28$. The SNRs of the real and the imaginary parts of the reconstruction are indicated in the subfigures.

Download Full Size | PDF

As plotted in Fig. 5(a) and 5(b) the SNRs of the real and imaginary parts of the reconstructed complex RI as functions of iteration number show the superiority of iODT to ODT. The nRMS errors in the amplitude and phase of the sinograms as functions of the iteration number are shown in Figs. 5(c) and 5(d), respectively, and they approximately converge to the floors corresponding to the solver discrepancy between FDTD and WABPM. After 28 iterations, the amplitude and phase of the field diffracted by the phantom and by the reconstructed object for the illumination angle of 0° are displayed in Figs. 5(e) and 5(f), respectively.

 

Fig. 5. SNRs of the real (a) and imaginary (b) parts of the reconstructed complex RI as functions of the iteration number m. nRMS errors in the amplitude (c) and phase (d) of sinograms as functions of the iteration number m. The dashed lines in (c) and (d) indicate the error floors corresponding to the solver discrepancy between FDTD and WABPM. The amplitude (e) and phase (f) of the field diffracted by the phantom (blue) and by the reconstructed object (red) after 28 iterations for the illumination angle of 0°.

Download Full Size | PDF

In the second example, we create a Shepp-Logan phantom with a higher RI contrast, as shown in Fig. 6(a). The real part of RI has a contrast of $({1.427 - {n_\textrm{b}}} )/{n_\textrm{b}} \approx 7\%$, and the imaginary part is identical to that in the previous example. Now, as shown in Fig. 6(b), conventional ODT reconstruction severely underestimates both the real and imaginary parts of complex RI. Furthermore, the crosstalk artifacts become much more pronounced in the imaginary part. The iODT reconstruction converges after $p = 27$ iterations. As shown in Fig. 6(c), $\textrm{Re}[{{n^{({27} )}}} ]$ is a good estimate of $\textrm{Re}[{{n_\textrm{t}}} ]$ but $\textrm{Im}[{{n^{({27} )}}} ]$ still has severe crosstalk artifacts. To suppress such artifacts, the proposed ES method is enabled for iteration numbers $m\;>\;27$. iODT with ES converges again at $m = 40$, where $\textrm{Re}[{{n^{({40} )}}} ]$ is very similar to $\textrm{Re}[{{n^{({27} )}}} ]$, and the reconstruction quality of $\textrm{Im}[{{n^{({40} )}}} ]$ is improved as displayed in Fig. 6(d). Note that when implementing iODT-ES, the real part of the reconstructed object function will be reset to the real part of the convergent iODT reconstruction, as shown in Algorithm 2 or Fig. 3. Since the imaginary part of the object function is much smaller than the real part, the real parts of reconstructed RI using iODT and iODT-ES, as shown in the upper row of Fig. 6(c) and 6(d), are almost identical.

 

Fig. 6. Distributions of the real (upper row) and imaginary (lower row) parts of the complex RI of (a) Shepp-Logan phantom with an index contrast of 7%, (b) ODT reconstruction, (c) iODT reconstruction at iteration $m = 27$, and (d) iODT-ES reconstruction at iteration $m = 40$, with ES module enabled for $m\;>\;27$. The SNRs of the real and the imaginary parts of the reconstruction are indicated in the subfigures.

Download Full Size | PDF

 

Fig. 7. SNRs of the real (a) and imaginary (b) parts of the reconstructed complex RI as functions of the iteration number m. nRMS errors in the amplitude (c) and phase (d) of sinograms as functions of the iteration number m. The ES module is enabled at iteration $m\;>\;27$. The dashed lines in (c) and (d) indicate the error floors set by the discrepancy between FDTD and WABPM simulations. The amplitude (e) and phase (f) of the field diffracted by the phantom (blue) and by the reconstructed object (red) after 40 iterations for the illumination angle of 0°.

Download Full Size | PDF

The SNRs of the real and imaginary parts of the reconstructed complex RI as functions of iteration number, as plotted in Fig. 7(a) and 7(b) respectively, show that iODT-ES provides the most accurate reconstruction compared to ODT and iODT. The nRMS errors in the amplitude and phase of the sinograms as functions of the iteration number are shown in Figs. 7(c) and 7(d), respectively. The reason why the sinogram errors (especially the amplitude error) suddenly increase at the 28^th iteration is that the initial guess for the iODT with the ES module is a phase object only, i.e., $\textrm{Re}[{{f^{({27} )}}} ]$, which produces non-absorptive diffracted fields. This is the same reason why the SNR in the imaginary part as shown in Fig. 7(b) suddenly increases at the 28^th iteration. The error curves in Figs. 7(c) and 7(d) approximately converge to the levels set by the solver discrepancy between FDTD and WABPM. After 40 iterations, the amplitude and phase of the field diffracted by the phantom and by the reconstructed object for the illumination angle of 0° are displayed in Figs. 7(e) and 7(f), respectively.

3.2 Phantom with photonic-crystal structure

To further assess the efficacy of iODT-ES, we use a phantom of a photonic-crystal fiber with 91 cores each of diameter 2.3 µm separated by 4.42 µm as shown in Fig. 8(a). The RI of the background medium is ${n_\textrm{b}} = 1.454$, which is perfectly matched with the RI of the cladding. The wavelength of the probe wave is $\lambda_0 = 532$ nm. The wavelength in the background medium is $\lambda_\textrm{b}{ = \lambda_0}/{n_\textrm{b}} = 365.9$ nm. The real part of the RI has a contrast $({1.469 - {n_\textrm{b}}} )/{n_\textrm{b}} \approx 1\%,$ and the imaginary part introduces an average energy loss of about 6% for the diffracted fields. The reconstruction dimension is ${170\lambda_\textrm{b}} \times {170\lambda_\textrm{b}}$ with a grid size ${\Delta}x = {\Delta}y{ = \lambda_\textrm{b}}/5$. Although each core has a relatively low contrast, an array of cores can still produce multiple-scattering effects. The object with the complex RI reconstructed by conventional ODT is shown in Fig. 8(b), in which the real part has incorrect values and the imaginary part contains severe artifacts. The iODT reconstruction (without the ES module) converges at the 35^th iteration as plotted in Fig. 8(c). Although the artifacts in the imaginary part become smaller compared with the ODT result (especially in the background region), it still has incorrect shapes and values for the lossy cores. To correct such artifacts, the proposed ES module is enabled for iterations $m\;>\;35$ and converges again at $m = 47$. As demonstrated in Fig. 8(d), iODT-ES subtracts out the error caused by the solver discrepancy, and therefore provides a more accurate reconstruction of the imaginary part. The data processing time is about one minute per iteration.

 

Fig. 8. Distributions of the real (upper row) and imaginary (lower row) parts of the complex RI of (a) phantom with photonic-crystal structure with an index contrast of 1%, (b) ODT reconstruction, (c) iODT reconstruction at iteration $m = 35$, and (d) iODT-ES reconstruction at iteration $m = 47$ with the ES module enabled for $m\;>\;35$. The SNRs of the real and the imaginary parts of the reconstruction are indicated in the subfigures.

Download Full Size | PDF

The SNRs of the real and imaginary parts of the reconstructed complex RI as functions of iteration number, as shown in Fig. 9(a) and 9(b) respectively, show that the reconstruction using iODT-ES is the most accurate compared to ODT and iODT. The nRMS errors in the amplitude and phase of the sinograms as functions of the iteration number are plotted in Figs. 9(c) and 9(d), respectively. As expected, the errors (especially in the amplitude) increase once the ES module is turned on at the 36^th iteration. The error curves in Fig. 9(c) and 9(d) do not approach the levels set by the solver discrepancy between FDTD and WABPM, partly because of the structural complexity of the object. Comparisons of the fields diffracted by the phantom and by the reconstructed object at the 47^th iteration for the illumination angle of 0° are shown in Figs. 9(e) and 9(f).

 

Fig. 9. SNRs of the real (a) and imaginary (b) parts of the reconstructed complex RI as functions of the iteration number m. nRMS errors in the amplitude (c) and phase (d) of sinograms as functions of the iteration number m. The ES module is enabled at iteration $m\;>\;35$. The dashed lines in (c) and (d) indicate the error floors set by the discrepancy between FDTD and WABPM. The amplitude (e) and phase (f) of the fields diffracted by the phantom (blue) and by the reconstructed object (red) after 47 iterations for the illumination angle of 0°.

Download Full Size | PDF

3.3 Noise tolerance

To test the efficacy of iODT-ES against noises, a Gaussian white noise (see Appendix) with a signal-to-noise ratio (SNR) of 20 dB is added to the simulated true fields diffracted by the photonic-crystal phantom shown in Fig. 8(a) and reproduced in Fig. 10(a). iODT reconstruction converges after 27 iterations. As shown in Fig. 10(b), the existence of noise does not degrade the reconstruction of the real part; however, there exist more severe artifacts in the imaginary part, from which the structure information can barely be identified. The ES module is then enabled after the 27^th iteration, and the reconstruction converges again at the 40^th iteration as displayed in Fig. 10(c). The reconstruction quality of the imaginary part improves in terms of artifact suppression, shape, and value.

 

Fig. 10. Distributions of the real (upper row) and imaginary (lower row) parts of the complex RI of (a) phantom with photonic-crystal structure that are identical to Fig. 8(a), (b) iODT reconstruction at iteration $m = 27$, and (c) iODT-ES reconstruction at iteration $m = 40$ with the ES module enabled for $m\;>\;27$. The SNRs of the real and the imaginary parts of the reconstruction are indicated in the subfigures.

Download Full Size | PDF

We tested the efficacy of ES module against different noise levels with SNRs ranging from 17 dB to 30 dB, along with the noise-free one as a reference. The nRMS errors in amplitude and phase increase as the noise level grows, as shown in Fig. 11. The phase error changes more slowly than the amplitude error, implying that the real-part reconstruction is less sensitive to noise or artifacts. The insets in Fig. 11 compare reconstructions of the imaginary part using iODT and iODT-ES for different noise levels, and demonstrate that the proposed ES module remains effective for SNR above 20 dB.

 

Fig. 11. nRMS errors in the amplitude and phase as functions of SNR. Insets are the imaginary part reconstructions using iODT and iODT-ES for different SNRs. The SNRs of the real and the imaginary parts of the reconstruction are indicated in the subfigures.

Download Full Size | PDF

4. Conclusion

Optical imaging of phase objects beyond the weakly-scattering regime is challenging, and the challenge is compounded by the presence of loss. While iterative optical diffraction tomography (iODT) enables successful reconstruction of non-lossy objects of higher contrast, compared to ODT, it is prone to artifacts in the presence of loss, much like ODT. In this study, we observed that for samples of high contrast, or complicated structure, inaccuracies in the forward-backwards solver (e.g., WABPM) produce field errors in the simulation volume, which are manifested in the filtered backpropagation algorithm as apparent crosstalk between the real and imaginary parts of the reconstruction. To help suppress these artifacts in samples whose object function is predominantly real, we assume that the magnitude of the error in the reconstructed object function, after iODT converges, is much smaller than the magnitude of the true object function. Under this premise, we may assume that the error in the real part of the object function is approximately negligible, while the error in the imaginary part (e.g., real-to-imaginary crosstalk) is not. We have exploited this differential property of the crosstalk effect by developing an error subtraction (ES) technique and designing an add-on module that reduces the artifacts of iODT. We have numerically demonstrated that iODT-ES is capable of simultaneously reconstructing complex refractive-index objects exhibiting multiple scattering. The new method is efficient and more accurate than iODT, and is robust in the presence of noise for SNR above 20 dB.