1. Introduction

Optical diffraction tomography (ODT) [1] is a mature phase-imaging technique used to reconstruct the 3D distribution of the complex refractive index. Because it is label-free, ODT is used in applications ranging from live cellular imaging [2] to refractive index (RI) profiling of optical fibers [3,4]. In conventional ODT, the sample to be imaged is illuminated with continuous-wave (CW) plane waves from multiple angles, and the respective diffracted optical fields are measured interferometrically. An inversion algorithm is then used to reconstruct the sample’s complex RI distribution.

An iterative Rytov-based ODT (iODT) algorithm [4] has recently been used for reconstruction of highly scattering objects. This algorithm works by forward propagating the known input field and backward propagating the measured output fields through an estimate of the sample’s RI distribution, for each illumination angle. The differential Rytov phase between these two fields is then used in a filtered backpropagation framework to calculate the error in the estimate of the sample. This error is then subtracted from the estimate to obtain a more accurate update of the RI distribution, and the process is repeated recursively.

One of the limitations of the iODT algorithm, under conditions of high-contrast RI distributions, is that the phase of the forward and backward propagated fields may become ill-defined in areas of the scattering volume where the fields contain small amplitudes. When this occurs, the Rytov phases used to reconstruct the RI become contaminated with phase-vortices that introduce challenges to correct phase unwrapping resulting in reconstruction artifacts. We have found that regularized optimization-based approaches that minimize a cost function based on the complex fields, or their respective complex phases [5–8], can be particularly effective at achieving accurate reconstructions for such objects because the iterative use of regularization (e.g., total-variation) alleviates challenges associated with poorly behaved signal, or missing information. Unfortunately, optimization algorithms that are based on a field fidelity criterion typically require a sufficiently accurate initial estimate to converge to a correct solution [8,9]. Furthermore, gradient-descent approaches may take several dozens of iterations before sufficiently converging to a solution. Unlike convex optimization approaches, however, iODT seeks a new “best estimate” of the sample’s RI distribution, rather than a gradual “descent” to the correct solution. This perturbative approach gives iODT an advantage in computational efficiency over optimization-based solutions.

In this paper, we introduce and validate a new strategy that combines iODT’s efficient framework with the benefits of regularized optimization techniques. This is achieved by applying total-variation (TV) regularization [5,10–12] in each iteration of the standard iODT algorithm. We call this method regularized iODT (R-iODT). We have validated this technique using simulated data and an experimental test and concluded that a substantial improvement in the signal-to-noise ratio (SNR) of reconstructions was obtained. Application of TV regularization at each iteration is significantly better than its application after the termination of the iteration.

2. R-iODT validation

2.1 iODT algorithm

The iODT algorithm [4], summarized in Alg. 1, uses a Rytov-like approach to iteratively reconstruct the RI distribution by using the complex-valued phase of diffracted fields. The “true,” or “measured,” diffracted fields are recorded holographically for each illumination angle in experimental tests [2,13,14], or numerically simulated by propagating the known input field through phantoms using a numerical solver such as the beam propagation method, or finite-difference time-domain (FDTD) method [15,16]. The object function, $f({x,y} )$, described in Alg. 1 is related to the RI by $f({x,y} )\equiv k_o^2[{{n^2}({x,y} )- n_b^2} ],$ where k_o is the wavenumber 2π/λ_o, and n_b is the background RI. As shown in Fig. 1, propagation is performed in local coordinates $({\xi ,\eta } )$, which relate to the global coordinates $({x,y} )$ by the transformation $x = \xi \cos \theta - \eta \sin \theta ,\; \; y = \xi \sin \theta + \eta \cos \theta .$ The true fields are measured along the local coordinate axis η, while propagation for each illumination angle θ is along the ξ axis.

 

Fig. 1. Schematic of an ODT experiment. Plane wave ${k_b}$ illuminates the sample at θ and travels along ξ direction. The scattered field is measured on the screen at $\xi = d$ along the η direction. The illumination angle is changed and the process is repeated.

Download Full Size | PDF

Using an initial guess of the sample f ⁽⁰⁾, the iODT algorithm takes a known input field u₀ for each illumination angle θ and forward propagates it through the sample to create $u_{\textrm{fwd}}^{(m )}({\xi ,\eta ,\theta } ).$ Likewise, the respective measured field for each illumination angle is backpropagated through the estimated RI to generate a field $u_{\textrm{bwd}}^{(m )}({\xi ,\eta ,\theta } )$. The forward and backward propagation is implemented using solvers, denoted by $\cal{S}$ and $\cal{S}$^-1, respectively. In this paper, we elect to use a “cold” initialization (${f^{(0 )}} = 0$) of the object’s distribution. For each iteration m, a differential “Rytov” phase, Δϕ_R is calculated and used to reconstruct the estimated error in the estimated distribution of the object $\mathrm{\Delta }{f^{(m )}},$ which is subtracted from the current estimate of the object ${f^{(m )}}$ to form a new estimate for the next iteration. The framework for 3-dimensional objects is straightforward but requires an updated Green’s function to perform field propagation and inversion. For computational brevity, we elect to perform reconstructions on two-dimensional objects, or 2-D cross-sections of extended objects (e.g., optical fibers).

Due to the multivalued nature of the algorithm used to calculate the Rytov phase, it is necessary to apply a phase-unwrapper is necessary on the imaginary part of Δϕ_R. For this study, we elect to use an L-2 norm phase-unwrapper [17] to perform this task. For objects that produce large amounts of scattering, the phase may contain local features vortices, where the value of the phase is ill-defined. Such vortices may introduce artifacts and even inhibit the process of unwrapping, thereby leading to reconstructive error.

At the end of each iODT iteration, constraints, denoted by the $\cal{C}$ operator, may be applied to the current estimate of the object (e.g., non-negativity). In this paper, we expand the iODT constraint operator to include TV-based regularization (e.g., TV), which is applied to the current estimate of the object at the end of each iteration. By doing so iteratively, the iODT algorithm can suppress artifacts caused by poor data quality, poorly behaved phase features due to phase unwrapping, or the presence of phase vortices.

This process is repeated until a stopping criterion is satisfied. For cross-comparison purposes, we elect to run the iODT (and R-iODT) algorithm(s) for a set number of iterations, m_max, which is specified for each experiment. The stopping condition, described in Alg. 1, is based on the convergence of the normalized root-mean-square (nRMS) errors in amplitude and phase of the sinogram:

2.2 TV algorithm

TV regularization is a popular regularization choice for convex optimization problems in tomographic phase imaging. One TV regularization approach, developed by Beck and Teboulle [12], and used in [5], employs a “dual approach” method to apply TV regularization on a reconstructed RI distribution. The modular nature of this approach is particularly useful for iODT since it accepts an input of an RI distribution, along with hyperparameters and constraints to return a new RI image with TV enhancement applied. We elect to adopt this modular TV algorithm so that both constraints on the RI distribution and regularization can be included in the $\cal{C}$ operator in Alg. 1 at the end of every iteration. This modular dual approach — specifically for isotropic TV in our study — is shown in Alg. 2 (from Appendix B in [5]), in which D is the discrete gradient operator applied to the input image (or datacube), proj_χ is a projection used to apply constraints, such that the values of the vectorized input RI distribution x are truncated to lie between n_min and n_max, ${[\mathbf{g} ]_n} \in {\mathrm{\cal{R}}^{3 \times \boldsymbol{N}}}$ is the gradient vector field of the discretized image $\boldsymbol{x} \in {\mathrm{\cal{R}}^{\boldsymbol{N}}}$ at the pixel position $n \in [{1, \ldots ,N} ],\boldsymbol{\; }$ and proj$_{\cal{G}}$ is the projection used in the case of isotropic TV:

3. R-iODT validation

In this section, we validate R-iODT through both simulations and experiments. To numerically validate R-iODT, reconstructions are performed on various phantoms that either exhibit complicated structure or large optical path-length difference (OPD). Data for the “true” sinogram fields is obtained using a FDTD solver with a grid spacing of λ_o/20, over a 360° span of illumination angles θ, in 5° increments. The illumination wavelength is λ_o= 1 µm. For each simulation, a reconstruction was performed using the standard iODT algorithm, with real-valued constraints on the RI distribution, and another using TV regularization, with the same constraints.

For samples whose RI distribution is known (i.e simulations in Sec. 3.1), the signal-to-noise-ratio

3.1 Validation through simulation

In the first example, Phantom 1, a simple distribution with high index contrast consisting of three disks, shown in Fig. 2(a), was used. Each disk has a 4.5 µm diameter, and RI of 1.348. The background RI ${n_b}$ is 1.518. For this simulation, phase unwrapping was turned off after 10 iterations. As shown in Fig. 2(b), while the phase unwrapping allowed iODT to obtain a maximum reconstructed SNR of 16.5 dB upon turning off phase unwrapping at the $m = 11$ iteration, phase vortices introduce artifacts at each subsequent iteration, which result in an overall degradation in SNR to 12.08 dB at iteration 50 (Fig. 2(b)). We note that ${\epsilon _A}$ diverges after 28 iterations for iODT, which can be directly attributed to the aforementioned vortex-introduced artifacts introduced at every iteration. The R-iODT reconstruction, however, shows a better SNR of 17 dB upon turning off phase unwrapping, and finishes with an SNR of 23.2 dB after 50 iterations. A maximum SNR of 22.9 dB is obtained using R-iODT, as shown in Fig. 2(b) at iteration $m = $ 38 and onwards, suggesting that iterative artifacts caused by phase vortices are suppressed, but not eliminated, by TV regularization, allowing the R-iODT algorithm to produce a more convergent and accurate solution. For objects containing well-defined regions of various RIs, TV regularization is keenly effective, as it preserves edges while suppressing total-variation in the background caused by artifacts.

 

Fig. 2. (a) Phantom 1 RI distribution, (b) SNR of the iODT and R-iODT reconstructions over the first 50 iterations. (c and d) nRMS errors in amplitude and phase of the sinogram ${\epsilon _A}$ and ${\epsilon _\phi }$ for iODT and R-iODT, over the first 50 iterations. Red arrows mark the iteration where the stopping criterion is satisfied for R-iODT (iODT diverges and so no iteration satisfies stop condition). (e and f) Respective iODT and R-iODT reconstructions at the m = 50 iteration.

Download Full Size | PDF

For Phantom 2, a more complex, walled, cell-like structure, shown in Fig. 3(a), was used with the background RI n_b = 1.518, n₁ = 1.378, n₂= 1.358, n₃ = 1.548, and n₄ = 1.568. The SNR of the reconstruction, shown in Fig. 3(b), shows that the R-iODT algorithm climbs to a progressively higher SNR reconstruction, while the iODT reconstruction obtains a maximum SNR of 14 dB at iteration 16, after which artifacts cause a steady degradation in reconstructive accuracy.

 

Fig. 3. (a) Phantom 2 RI distribution, (b) SNR of the iODT and R-iODT reconstructions over the first 50 iterations. (c and d) nRMS errors in amplitude and phase of the sinogram ${\epsilon _A}$ and ${\epsilon _\phi }$ for iODT and R-iODT, over the first 50 iterations. Blue and red arrows mark the iteration where the stopping criterion is satisfied. (e and f) Respective reconstructions after stopping criterion is satisfied (iteration $m = \; 43$ for iODT and $m = \; 49$ for R-iODT, respectively). (g and h) Respective iODT and R-iODT reconstructions after m = 50 iterations.

Download Full Size | PDF

The nRMS errors of the output fields’ amplitude and phase are plotted in Fig. 3(c) and 3(d), respectively, and show that both algorithms yield reconstructions that reduce error in output fields, even though the iODT reconstruction contains a larger SNR error. This suggests that regularization can guide the algorithm away from possible solutions that presumably contain higher total variance than the correct distribution. Using $Q = 10$ for the stopping criterion outlined in 2.1, based on the ${\epsilon _A}$ and ${\epsilon _{ph}}$ plots in Fig. 3(c) and 3(d), respectively, the iODT algorithm would have satisfied the stopping criterion in Alg. 1 at the $m = 43$ iteration, and the R-iODT algorithm at the $m = 49$ iteration, with corresponding reconstructions shown in Fig. 4(e) and 4(f). The iODT and R-iODT reconstructions after the m = 50 iteration are shown in Fig. 3 g and 3 h.

 

Fig. 4. (a) Phantom 3 RI distribution, (b) SNR of the iODT and R-iODT reconstructions over the first 50 iterations. (c and d) nRMS errors in amplitude and phase of the sinogram for both algorithms. (e and f) Respective reconstructions after stopping condition is satisfied (iteration $m = \; 27$ for iODT and $m = \; 48$ for R-iODT). (g and h) Respective iODT and R-iODT reconstructions after $m = \; 50$ iterations.

Download Full Size | PDF

The plots in Fig. 3(b)-(d) exhibit a non-smooth feature at the m = 10 iteration, marking the iteration where phase-unwrapping is turned off. Unlike iODT, however, a noticeable drop in the ${\epsilon _A}$ and ${\epsilon _{ph}}$ nRMS errors for R-iODT near iteration 30, suggesting that the algorithm was able to “lock on” to a solution that better matched both the output fields and the correct RI distribution. The standard iODT algorithm, however, shows a decline in both nRMS error and SNR, suggesting that the algorithm was trapped in a solution whose output fields that match the measured ones, but contains the incorrect RI distribution.

As demonstrated in Figs. 2(c) and 3(e) for Phantoms 1 and 2, the iODT algorithm is sensitive to phase-vortex contamination of the Rytov phase. These vortices can appear in objects with large OPDs, as well as objects of complicated structure, whose features contain large RI contrast. To demonstrate the effectiveness of R-iODT for such objects, reconstructions are performed on Phantom 3, a 19-disk distribution shown in Fig. 4(a), where n₁= 1.38, n₂= 1.418, and n_b= 1.518. Although the overall OPD of the object is less than 2π, the complicated structure of many small, high-contrast features produces phase vortices that contaminate the Rytov phase. Although iODT and R-iODT have similar SNR upon turning off phase unwrapping at iteration 11, phase vortex artifact contamination limits the SNR of the iODT reconstruction to 4.7 dB. Like the previous objects, the contaminated phase introduces error at each iODT iteration, which begin to lower the SNR of the iODT reconstruction after iteration 22. For R-iODT, these artifacts are suppressed at each iteration, allowing the algorithm to reach an SNR of 14.3 dB. We note that although the convergence condition is satisfied for both algorithms, the SNR of the reconstructions continue to change even at the 50^th iteration, suggesting that the stopping criterion might also include a condition based on values in the RI distribution converging between consecutive iterations, in addition to the errors in Eqs. (1) and (2).

For the three simulations above, each iODT iteration takes on average of 35 seconds to complete. The R-iODT algorithm differs only in the use of a TV-enhancement step at the end of each iteration. The latter method, on average takes 42 seconds to complete. The code for the iODT algorithm, however, has been optimized for the CPU and unparallelized computation, and much shorter iteration times are possible.

To assess the feasibility of reducing the computational burden of R-iODT versus the standard iODT, we compare the SNR of R-iODT (TV regularization at each iteration) to a reconstruction where TV regularization is performed only once, after the final iODT reconstruction. The result, shown in Fig. 5(c), has an SNR of 4.8 dB, demonstrating that the R-iODT algorithm (Fig. 5(b)) reconstructs more accurately when TV regularization is applied at each iteration (or every few iterations) rather than after the final iODT reconstruction.

 

Fig. 5. Phantom 3 reconstructions (a) iODT, (b) R-iODT using TV at each iteration, and (c) iODT with TV regularization applied to final $({m = 50} )$ reconstruction.

Download Full Size | PDF

While the numerical validation of R-iODT, using TV regularization, has shown considerable improvement over the standard iODT algorithm, the choice of the type of regularization used for the type of sample must be carefully made. When reconstructing features that are of similar scale to the resolution limit of ODT (i.e., ${\sim} {\lambda _o}/2$), the TV enhancement may lower the overall SNR of small, reconstructed features if the algorithm is unable to sharpen such features while preserving their correct values. For such samples, a more sophisticated regularization approach may be considered.

3.2 Experimental validation

We have applied the R-iODT algorithm to holographically measured ODT data taken on a 19-core step-index multicore fiber using an Intrafiber IFA-100 optical fiber profiler at λ₀ = 630 nm. The raw dataset used is the same as in [4]. Reconstructions using the standard iODT without TV regularization, applying TV regularization once after the final standard iODT iteration, and using R-iODT are shown in Fig. 6. A cleaner reconstruction in the central region of the fiber is observed when iteratively applying regularization in the R-iODT algorithm. The pixel size for the reconstruction is Δx = 0.184 µm, and the diffracted fields are measured on the object over 36 angles, ranging from 0° to 175°, in intervals of 5°.

 

Fig. 6. Reconstructions of multicore optical fiber (a) SEM of fiber, (b) iODT, (c) iODT with TV regularization applied to final iteration $({m = 50} )$, and (d) R-iODT using TV at each iteration.

Download Full Size | PDF

4. Conclusion

Iterative optical diffraction tomography (iODT) is a perturbative approach that offers an advantage in computational efficiency over general optimization-based solutions, for large contrast objects [9]. We have developed a new version of iODT, named R-iODT, in which total-variation (TV) regularization is applied to the RI estimate at each iteration. Numerical and experimental tests demonstrated significant improvement in SNR of the reconstructed RI and confirmed reduction of the error between the reconstructed fields and the true fields. Regularization at each iteration tends to suppress the effects of artifacts that may occur in the case of highly-scattering samples. The R-iODT algorithm accurately and efficiently yields reconstruction of high-contrast samples without heavily relying on the careful selection of hyper-parameters or run-over large numbers of iterations that are encountered in general optimization approaches.

The application of TV regularization at each iteration was found to be significantly superior to one time application at the final iteration of a standard iODT reconstruction. Regularization seems to guide the iterative process toward better reconstruction. On average, the additional computational time per iteration increased by ∼20% compared to the standard iODT algorithm. The cost to computational time for 3-dimensional RI distributions is expected to produce proportionally similar increase in computational time for 2-D objects.

In this paper, we have limited our tests of R-iODT to TV regularization. Other regularization techniques can also be used and might offer better improvements in SNR, and possibly resolution enhancement. Other applications for regularization, such as alleviating artifacts caused by experimental limitations (e.g., missing angles problem) would also be valuable to explore in future work.