1. Introduction

Optical diffraction tomography (ODT) was introduced in [1] by E. Wolf in the late ’60s. It is a microscopic technique that retrieves the distribution of refractive indices in biological samples out of holographic measurements of the scattered complex field produced when the sample is illuminated by an incident wave. This method is of particular interest in biology because, contrarily to fluorescence imaging, it does not require any staining of the sample [2]. It proceeds by solving an inverse scattering problem, where the scattering phenomenon is governed by the wave equation. There is a vast literature on inversion methods going from linearized models (Born, Rytov) [1, 3] to nonlinear ones [4–7]. It is worth noting that the scattering model, along with its associated inverse problem, is generic and not limited to optical diffraction tomography. In particular, it is encountered in many other fields such as acoustics, microwave imaging, or radar applications [8].

1.1. From the wave equation to the Lippmann-Schwinger integral equation

Let us consider an unknown object of refractive index n(x) lying in the region Ω ⊆ ℝ^D (D ∈ {2, 3}) and being immersed in a medium of refractive index n_b, as depicted in Fig. 1. This sample is illuminated by the incident plane wave

 

Fig. 1 Optical diffraction tomography. A sample of refractive index n(x) is immersed in a medium of index n_b and illuminated by an incident plane wave (wave vector k). The interaction of the wave with the object produces forward and backward scattered waves. The forward scattered wave is recorded in the detector plane. Optionally, a second detector plane may record the backward scattered wave (see Section 5).

Download Full Size | PDF

1.2. Inverse ODT problem: prior work

Let the object be illuminated by a series of incident fields $u_p^{\text{in}}(\mathbf{x})$, p ∈ [1 … P]. Records of the resulting total fields u_p (x) at positions x_m (m ∈ [1 … M]) in the detector plane Γ are denoted y_p ∈ ℂ^M (see Fig. 1). The objective is then to retrieve the scattering potential function f (x) (i.e., the refractive index n(x)) from the data y_p. Pioneering methods were using linear approximations of the model. For instance, assuming that the scattering field is weak compared to the incident one (i.e., u^sc ≪ uⁱⁿ), one can interpret the phase of the transmitted wave as the Radon transform of the refractive index and then reconstruct f using the filtered-back-projection algorithm [10, 11]. This method ignores the effect of diffraction. The first Born approximation [1] has then been proposed as a refined model. Its validity is however limited to thin samples with weak variations of their refractive index (RI) [12]. A more accurate linearization, less sensitive to the thickness of the sample but still limited to weak RI contrasts, is given by the Rytov approximation [3, 13]. It is derived by assuming that the total field has the form u(x) = uⁱⁿ(x)e^ϕ(x), where ϕ(x) is a complex phase function. Both Born and Rytov approximations have been originally used to derive direct inversion methods. They were later used within regularized variational approaches to improve the quality of reconstructed images [14,15].

Inversion methods that use a nonlinear model have been shown to significantly improve the accuracy of reconstruction. These include the conjugate-gradient method (CGM) [16, 17], the contrast source-inversion method (CSI) [18], the beam-propagation method (BPM) [19], the recursive Born approximation [5], or the hybrid method proposed in [4]. Although still approximate (for instance, they do not properly take reflections into account), they more closely adhere to the model of the physical phenomenon than the linear models, at the price of a higher computational cost. We refer the reader to [2] for additional details concerning existing approximations, regularizations, algorithms, and comparisons.

To address applications with thick samples and large RI contrasts, a better solution is to rely on the exact Lippmann-Schwinger model which accounts for mutiple scattering and reflections. Such an approach has been recently proposed in [6,7] (SEAGLE algorithm). There, the authors tackle the problem from a variational perspective. They minimize a nonconvex objective using the well known fast iterative shrinkage-thresholding algorithm (FISTA) [20]. Their main contribution is to compute the forward model (which itself requires the inversion of an operator) using Nesterov’s accelerated gradient-descent (NAGD) method [21] and, more interestingly, to explicitly compute the gradient of the quadratic data-fidelity term as an error-backpropagation of the forward algorithm. However, the bottleneck of their method is its high memory consumption. Indeed, the error-backpropagation strategy requires one to store all the iterates produced during the computation of the iterative forward model. This can be limiting for large 3D volumes.

1.3. Contributions

To improve the computational efficiency of solvers such as SEAGLE, we provide an explicit expression for the Jacobian of the nonlinear Lippmann-Schwinger operator. This results in an efficient method to compute the gradient of the data-fidelity term and avoid recoursing to the memory-consuming error-backpropagation strategy. Another advantage is that the computation of the forward model and of the gradient are now decoupled. They can thus be solved using any numerical scheme. Then, considering simulated data, we show that the proposed method results in a significant reduction of both computational time and memory requirements with respect to SEAGLE, at no loss in quality.

In Section 2.1, we formulate the discrete forward model proposed in [7]. Then, the common approach used to solve the inverse problem subject to sparsity constraints is presented in Section 2.2. There, we highlight our main innovation with respect to SEAGLE, which is a new computation of the gradient of the data-fidelity term. It relies on the derivation of the Jacobian of the forward model, which is given by Proposition 3.1. Finally, Sections 4 and 5 are dedicated to numerical comparisons.

2. Solving the inverse problem

2.1. Formulation of the forward model

In this section, we review the formulation of the forward model that was proposed by Liu et al. in [7]. Let the region of interest Ω be divided into N ∈ ℕ “pixels”. Then, over Ω, we define the discrete version of Eq. (7) as

2.2. Common optimization strategy

Following the classical variational approach, the estimation of f ∈ ℝ^N from the measurements {y_p ∈ ℂ^M}_{p∈[1 … P]} is formulated as the optimization problem

Algorithm 1. Accelerated forward-backward splitting.

View Table | View all tables in this article

2.2.1. Computation of ${\mathbf{J}}_{h_p}^H(\mathbf{f})$

The computation of ${\mathbf{J}}_{h_p}^H(\mathbf{f})$, required at line 5 of Algorithm 1, is challenging. The existence of a closed-form solution is made unlikely by the fact that the forward model in Eq. (8) itself requires one to invert an operator. We distinguish two distinct strategies.

2.2.2. Computation of ${\text{prox}}_{\gamma \lambda ℛ}$

Numerous methods have been proposed to compute the proximity operator of $ℛ$, [25–27]. In SEAGLE, Liu et al. use the algorithm proposed by Beck and Teboulle [25]. Here, we compute it using the popular alternating-direction method of multipliers (ADMM) [28], which is well suited to the minimization of the sum of three convex functions. Moreover, it provides a high modularity for spatial regularization since one can easily change from one regularizer (e.g., TV) to another (e.g., Hessian Shatten-norm [29]). Details about the computation of $ℛ$ are provided in Appendix A.

2.2.3. Speedup strategies

The cost of evaluating the forward model with Eq. (9) and the gradient $\mathbf{\nabla }\mathcal{D}$ is proportional to the number P of illuminations ${\mathbf{u}}_p^{in}$. However, these computations can easily be parallelized by performing the computation for each illumination (or each element of the sum in Eq. (15)) on a separate thread. Moreover, in the spirit of the stochastic gradient-descent algorithm [30], we approximate $\mathbf{\nabla }\mathcal{D}$ as

3. Efficient computation of the gradient $\mathbf{\nabla }\mathcal{D}$

The error-backpropagation strategy used in SEAGLE to compute ${\mathbf{J}}_{h_p}^H(\mathbf{f})$ implies that one must store all the forward iterates. This consumes memory resources and compromises the deployment of the method for large 3D data. Instead, Proposition 3.1 reveals that its computation requires one to invert the operator (I − diag(f)G^H). This operator has the same form (and size) that the operator we invert within the forward computation in Eq. (9) and both can be computed in a similar way, using an iterative algorithm. Moreover, it allows us to decouple the forward and gradient computation in Algorithm 1, which has the two following advantages:

Proposition 3.1. The Jacobian matrix of the function h_p in Eq. (17) is given by

Then, from Eq. (8), we get that

4. Algorithm analysis

4.1. Memory requirement

In this section, we elaborate on the memory consumption of the proposed method in comparison with SEAGLE. First, let us state that gradient based methods, such as NAGD or CG, have similar memory requirements. It corresponds roughly to three times the size of the optimization variable which is the part that is common to both algorithms. The additional memory requirement that is specific to SEAGLE relies only on the storage of the NAGD iterates during the forward computation. Suppose that K_NAGD ∈ ℕ iterations are necessary to compute the forward model with Eq. (9) and that the region Ω is sampled over N ∈ ℕ pixels (voxels, in 3D). Since the total field u_p (f) computed by NAGD is complex-valued, each pixel is represented with 16 bytes (double precision for accurate computations). Hence, the difference of memory consumption between SEAGLE and our method is

For illustration, we give in Fig. 2 the evolution of Δ_Mem as a function of N for different values of K_NAGD and N_Threads. One can see with the vertical dashed lines that, for 3D volumes, the memory used by SEAGLE quickly reaches several tens of Megabytes, even for small volumes (e.g., 128 × 128 × 128), to hundreds of Gigabytes for the larger volumes that are typical of microscopy (e.g., 512 × 512 × 256). This shows the limitation of SEAGLE for 3D reconstruction in the presence of a shortage of memory resources and reinforces the interest of the proposed alternative.

 

Fig. 2 Predicted evolution of Δ_Mem as function of the number N of points for two values of K_NAGD and N_Threads. The vertical dashed lines give examples of 2D and 3D volumes for a range of values of N. Finally, the three crosses correspond to values of Δ_Mem measured experimentally.

Download Full Size | PDF

4.2. Conjugate gradient vs. Nesterov accelerated gradient descent for Eq. (9)

Due to Proposition 3.1, we can compute both Eq. (9) and $J_{h_p}^H(\mathbf{f})$ using any state-of-the-art quadratic optimization algorithm. This contrasts with SEAGLE, where one must derive the error-backpropagation rule from the forward algorithm, which may limit its choice. We now provide numerical evidence that GC is more efficient than NAGD for solving Eq. (9). To this end, we consider a circular object (bead) of radius r_bead and refractive index n_bead immersed into water (n_b = 1.333), as presented in Fig. 3 (top-left). In such a situation, an analytic expression of the total field is provided by the Mie theory [31,32]. Hence, at each iteration k, we compute the relative error ε_k of the current estimate u^k to the Mie solution u_Mie as

 

Fig. 3 Forward-model solution for a bead with radius 3λ and a contrast of 1 using CG (bottom-left) and NAGD (bottom-right), as well as the Mie solution (top-right). The setting used for this experiment is presented in the top-left panel. The colormap is the same for each figure.

Download Full Size | PDF

 

Fig. 4 Evolution of the number of iterations $k_{{\epsilon }_0}$ needed to let the relative error Eq. (28) fall below ε₀ = 10⁻² as function of bead radius (left) and bead contrast (right).

Download Full Size | PDF

Finally, the solution obtained with the two algorithms for r_bead = 3λ and a contrast of 1 are shown in Fig. 3. The analytic Mie solution is also provided for comparison. From these figures, one can appreciate the high accuracy obtained by solving Eq. (9), as first demonstrated in [6, 7].

5. Numerical experiments

This section is devoted to numerical experiments that illustrate the two main advantages of the proposed method over SEAGLE, which consist in a reduced computational cost and a reduced memory consumption. The algorithms have been implemented using an inverse-problem library developed in our group [33] (GlobalBioIm: http://bigwww.epfl.ch/algorithms/globalbioim/). Hence, they share the implementation of the overall FISTA algorithm as well as inner procedures such as the computation of the proximity operator of ℛ (see Appendix A). The only difference between the two methods resides in the computation of the forward model in Eq. (9) and $J_{h_p}^H(\mathbf{f})$. For SEAGLE, this is performed using the NAGD algorithm and an error-backpropagation strategy. For our method, Eq. (9) and $J_{h_p}^H(\mathbf{f})$ are computed using the CG algorithm, in accordance with Proposition 3.1. Note that no parallelization is used. Reconstructions are performed with MATLAB 9.1 (The MathWorks Inc., Natick, MA, 2000) on a PowerEdge T430 Dell computer (Intel Xeon E5-2620 v3).

5.1. Simulated data

5.1.1. Simulation settings

The Shepp-Logan phantom of Fig. 5 has the contrast $\max (\left|\mathbf{f}\right|)/\left(k_0^2{n}_{\mathrm{b}}^2\right)=0.2.$ It is immersed into water (n_b = 1.333). The wavelength of the incident plane waves is λ = 406 nm. We consider thirty-one incident angles, from −60° to +60^°. The sources are placed at the bottom side of the sample, at a distance of 16.5λ from its center. Moreover, we consider two detectors placed on both top and bottom sides of the object, also at a distance of 16.5λ from its center. Hence, the overall region is a square of length 33λ per side. Data are simulated using a fine discretization of this region, with a (1024 × 1024) grid that leads to square pixels of surface (3.223 10⁻²λ)². We used a large number of CG iterations to get the accurate simulation mentioned in Section 4.2. Then, the measurements were extracted from the first and last rows of each total field associated to the incident fields. This lead to a total of (31 × 2 × 1024) measurements. Finally, we defined three ODT problems by downsizing (using averaging) the (31 × 2 × 1024) measurements to grids with size of (31 × 2 × 512), (31 × 2 × 384), and (31 × 2 × 256).

 

Fig. 5 Sheep-Logan phantom and refractive indices of the gray levels. The contrast is 20%.

Download Full Size | PDF

This setting corresponds to an ill-posed and highly scattering situation. Moreover, the detector length is only two times larger than the object, which results in a loss of information for large incident angles. This makes the resulting inverse problem challenging.

5.1.2. Algorithm parameters

For each simulated OTD problem, we considered a square region of interest Ω with sides half the sources–detector distance. That corresponds to images of size (256 × 256) with pixels of area (6.445 · 10⁻²λ)², (192 × 192) with pixels of area (8.839 · 10⁻²λ)², and (128 × 128) with pixels of area (1.289 · 10⁻¹λ)². The support of the phantom is·fully contained in Ω.

Then, to compute the gradient (stochastic-gradient strategy), we selected eight angles over the thirty-one that were available and changed this selection at each iteration (see Section 2.2.3).

The NAGD or CG forward algorithms are stopped either after hundred-twenty iterations or when the relative error between two iterates is below 10⁻⁴. Finally, two-hundred iterations of FISTA are performed with a descent step fixed empirically to γ = 5 · 10⁻³. We used the regularization parameter μ = 3.3 · 10⁻².

5.1.3. Metrics

We compared the two methods in terms of running time and memory consumption, as measured by the peak memory (maximum allocated memory) reached by each algorithm during execution. The outcome is reported in Table 1. Once again, due to the use of our inverse-problem library [33], the comparison of the two methods is fair because their implementations differ only by the forward algorithm and by the computation of $J_{h_p}^H(\mathbf{f})$. Moreover, CG and NAGD are implemented in the same fashion since they inherit the same optimization class of our inverse-problem library. Finally, we also provide the SNR of the reconstructed refractive index and observe that the computational gain comes at no cost in quality.

Table 1. Proposed method vs. SEAGLE [6, 7] in terms of running time and memory consumption. The reconstructed refractive-index maps are presented in Fig. 6.

View Table | View all tables in this article

5.1.4. Discussion

Our proposed alternative to SEAGLE allows us to reduce both time and memory. Moreover, we have measured the peak memory difference Δ_Mem between the two methods and superimposed it on the predictions of Fig. 2 where the adequacy between the theoretical curves and the measured points is remarkable. Hence, although our experiments are restricted to 2D data, where the gap between the two algorithms is moderate, the evolution of Δ_Mem for 3D data can be extrapolated from Fig. 2. This shows the relevance of our method when size increases.

The SNR values given in Table 1 as well as the reconstructions presented in Fig. 6 suggest that the two methods perform similarly in terms of quality. This is not surprising since the overall algorithm is the same, the differences residing merely in the computation of the forward model in Eq. (9) and the Jacobian $J_{h_p}(\mathbf{f})$. Moreover, one can observe that the quality of reconstruction decreases when the discretization grid becomes coarser. Indeed, the model is insufficiently accurate when the discretization is too poor. For instance, in the case of the (128 × 128) grid, one wavelength unit is discretized using eight pixels, which is clearly detrimental to the accuracy of the forward model.

 

Fig. 6 Reconstructions obtained by the proposed method and by SEAGLE for the (256 × 256) ODT problem with μ = 3.3 · 10⁻². The colormap is identical to that of Fig. 5. For comparison, we provide the TV-regularized Rytov reconstruction with μ = 3 · 10⁻³.

Download Full Size | PDF

Reconstructions for the (256×256) problem are presented in Fig. 6 for completeness. Besides the difficulty of the considered scenario, the two methods are able to retrieve most details of the object in comparison with the Rytov approximation. Artifacts are mainly due to the missing-cone problem and to the limited length of the detector. This corroborates the findings of [7].

5.2. Real data

We evaluated our method using the FoamDielExt target (TM polarisation) of the Institut Fresnel’s public database [34]. The data were collected for the two-dimensional inhomogeneous sample depicted in the left panel of Fig. 7. The permitivity of the ground truth was measured experimentally and is subject to uncertainties [34]. The object is fully contained in a square region of length 15 cm per side, which we discretize using a 256 × 256 grid. Sensors were distributed circularly around the object, at a distance of 1.67 m from its center, and with a step of 1°. Eight sources, uniformly distributed around the object, were sequentially activated. For each activated source, the sensors closer than 60° from the source were excluded. Thus, 241 detectors among the 360 available were used for each source. Frequencies from 2 to 10 GHz with a step of 1 GHz are available in the database but we used only the 3 GHz measurements (i.e., λ = 10 cm).

 

Fig. 7 Reconstructions (permittivity) obtained by the proposed method and by SEAGLE for the FoamDielExt target of the Institut Fresnel’s database [34] with μ = 1.6 · 10⁻². The SNR values (computed from the experimentally measured permittivity of the ground truth) are 25.13 dB (Ours) and 25.15 dB (SEAGLE) while the computing times are respectively of 6 min and 93 min.

Download Full Size | PDF

The NAGD or CG forward algorithms are stopped either after two-hundred iterations or when the relative error between two iterates is below 10⁻⁶. Hundred iterations of FISTA are performed with a descent step γ = 5 · 10⁻³. We used the regularization parameter μ = 1.6 · 10⁻².

In Fig. 7, we see that both methods provide good reconstructions that are essentially indistinguishable (see also SNR values provided in the caption of the figure). This corroborates the simulated numerical experiments of Section 5.1. The main point here is that, for this setting, the proposed method was 15 times faster than SEAGLE.

6. Conclusion

We have presented a refinement of the SEAGLE algorithm that was recently proposed in [6, 7] and that has shown unprecedented reconstructions for difficult configurations. However, the current limitation of SEAGLE is that its memory requirements increase excessively with the size of the problem, particularly in 3D. As an alternative, we have derived the explicit expression of the Jacobian matrix $J_{h_p}(\mathbf{f})$ of the nonlinear Lippman-Schiwnger model and shown that it can be computed in a direct analogy with the computation of the forward model. This approach allows us to drastically reduce the memory consumption and opens the door to 3D reconstruction using desktop computers. Moreover, the proposed method is quite flexible in the sense that it can cope with any iterative algorithm employed to compute either the forward model or $J_{h_p}^H(\mathbf{f})$. For instance, the conjugate-gradient algorithm proved its efficiency for this task. It allows a significant decrease of the computational time with respect to SEAGLE. Finally, these improvements in terms of speed and memory come at no loss in quality.

A. Proximity operator of ℛ

In this appendix, we describe how we compute the proximity operator of the regularization term ℛ in Eq. (14) using the ADMM algorithm [28]. The proximity operator is defined [35] as the solution of the optimization problem

(29)$${\text{prox}}_{\mu ℛ}(\mathbf{v})=\underset{f\in {ℝ}^N}{\arg \min }\left(\frac{1}{2}{\Vert \mathbf{f}-\mathbf{v}\Vert }_2^2+\mu {\Vert \mathbf{f}\Vert }_{\text{TV}}+i{⩾}_0(\mathbf{f})\right)$$

for μ > 0. Let us start by reformulating Eq. (29) as

(30)$$\begin{array}{l}{\text{prox}}_{\mu ℛ}(\mathbf{v})=\underset{f\in {ℝ}^N}{\arg \min }\left(\frac{1}{2}{\Vert \mathbf{f}-\mathbf{v}\Vert }_2^2+\mu {\Vert {\mathbf{q}}_1\Vert }_{2,1}+i{⩾}_0({\mathbf{q}}_2)\right),\\ s.t.{\mathbf{q}}_1=\mathbf{\nabla }\mathbf{f},\\ {\mathbf{q}}_2=\mathbf{f},\end{array}$$

which admits the augmented-Lagrangian form

(31)$$\begin{array}{l}ℒ(\mathbf{f},{\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{w}}_1,{\mathbf{w}}_2)=\frac{1}{2}{\Vert \mathbf{f}-\mathbf{v}\Vert }_2^2+\frac{{\rho }_1}{2}{\Vert \mathbf{\nabla }\mathbf{f}-{\mathbf{q}}_1+\frac{{\mathbf{w}}_1}{{\rho }_1}\Vert }_2^2\\ +\frac{{\rho }_2}{2}{\Vert \mathbf{f}-{\mathbf{q}}_2+\frac{{\mathbf{w}}_2}{{\rho }_2}\Vert }_2^2+\mu {\Vert {\mathbf{q}}_1\Vert }_{2,1}+i{⩾}_0({\mathbf{q}}_2),\end{array}$$

where ρ₁ and ρ₂ are positive scalars, and where w₁ ∈ ℝ^{N ×D} and w₂ ∈ ℝ^N are the Lagrangian multipliers. Then, one can minimize Eq. (31) using ADMM. The iterates are summarized in Algorithm 2.

Algorithm 2. ADMM for solving Eq. (29).

View Table | View all tables in this article

For the sake of completeness, we provide in Eq. (32) and Eq. (33) the expressions

(32)$$\forall \mathbf{q}\in {ℝ}^N,{\left[{\text{prox}}_{i{⩾}_0}(\mathbf{q})\right]}_n={({\mathbf{q}}_n)}_{+},$$

(33)$$\forall \mathbf{q}\in {ℝ}^{N\times D},{\left[{\text{prox}}_{\gamma {\Vert \cdot \Vert }_{2,1}}(\mathbf{q})\right]}_{n,d}={\mathbf{q}}_{n,d}{\left(1-\frac{\gamma }{{\Vert {\mathbf{q}}_{n,.}\Vert }_2}\right)}_{+},$$

of ${\text{prox}}_{i{⩾}_0}$ and ${\text{prox}}_{\gamma {\Vert \cdot \Vert }_{2,1}}$

(34)$${(x)}_{+}:=\max (x,0),x\in ℝ.$$

Funding

European Research Council (ERC) (692726)

Acknowledgments

This research was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme, Grant Agreement no 692726 “GlobalBioIm: Global integrative framework for computational bio-imaging.”

References and links

1. E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Commun. 1, 153–156 (1969). [CrossRef]  

2. D. Jin, R. Zhou, Z. Yaqoob, and P. T. So, “Tomographic phase microscopy: Principles and applications in bioimaging,” J. Opt. Soc. Am. B 34, B64–B77 (2017). [CrossRef]  

3. A. Devaney, “Inverse-scattering theory within the Rytov approximation,” Opt. Lett. 6, 374–376 (1981). [CrossRef]   [PubMed]  

4. E. Mudry, P. C. Chaumet, K. Belkebir, and A. Sentenac, “Electromagnetic wave imaging of three-dimensional targets using a hybrid iterative inversion method,” Inverse Probl. 28, 065007 (2012). [CrossRef]  

5. U. S. Kamilov, D. Liu, H. Mansour, and P. T. Boufounos, “A recursive Born approach to nonlinear inverse scattering,” IEEE Signal Process. Lett. 23, 1052–1056 (2016). [CrossRef]  

6. H.-Y. Liu, U. S. Kamilov, D. Liu, H. Mansour, and P. T. Boufounos, “Compressive imaging with iterative forward models,” in “IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP),” (IEEE, 2017), pp. 6025–6029.

7. H.-Y. Liu, D. Liu, H. Mansour, P. T. Boufounos, L. Waller, and U. S. Kamilov, “SEAGLE: Sparsity-driven image reconstruction under multiple scattering,” arXiv preprint arXiv:1705.04281 (2017).

8. D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, vol. 93 (Springer Science & Business Media, 2012).

9. J. A. Schmalz, G. Schmalz, T. E. Gureyev, and K. M. Pavlov, “On the derivation of the Green’s function for the Helmholtz equation using generalized functions,” Am. J. Phys. 78, 181–186 (2010). [CrossRef]  

10. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (SIAM, 2001). [CrossRef]  

11. W. Choi, C. Fang-Yen, K. Badizadegan, S. Oh, N. Lue, R. R. Dasari, and M. S. Feld, “Tomographic phase microscopy,” Nat. Methods 4, 717 (2007). [CrossRef]   [PubMed]  

12. B. Chen and J. J. Stamnes, “Validity of diffraction tomography based on the first Born and the first Rytov approximations,” Appl. Opt. 37, 2996–3006 (1998). [CrossRef]  

13. Y. Sung, W. Choi, C. Fang-Yen, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Optical diffraction tomography for high resolution live cell imaging,” Opt. Express 17, 266–277 (2009). [CrossRef]   [PubMed]  

14. Y. Sung and R. R. Dasari, “Deterministic regularization of three-dimensional optical diffraction tomography,” J. Opt. Soc. Am. A 28, 1554–1561 (2011). [CrossRef]  

15. J. Lim, K. Lee, K. H. Jin, S. Shin, S. Lee, Y. Park, and J. C. Ye, “Comparative study of iterative reconstruction algorithms for missing cone problems in optical diffraction tomography,” Opt. Express 23, 16933–16948 (2015). [CrossRef]   [PubMed]  

16. P. C. Chaumet and K. Belkebir, “Three-dimensional reconstruction from real data using a conjugate gradient-coupled dipole method,” Inverse Probl. 25, 024003 (2009). [CrossRef]  

17. K. Belkebir, P. C. Chaumet, and A. Sentenac, “Superresolution in total internal reflection tomography,” J. Opt. Soc. Am. A 22, 1889–1897 (2005). [CrossRef]  

18. A. Abubakar and P. M. van den Berg, “The contrast source inversion method for location and shape reconstructions,” Inverse Probl. 18, 495 (2002). [CrossRef]  

19. U. Kamilov, I. Papadopoulos, M. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Learning approach to optical tomography,” Optica 2, 517–522 (2015). [CrossRef]  

20. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2, 183–202 (2009). [CrossRef]  

21. Y. Nesterov, “A method of solving a convex programming problem with convergence rate O(1/k²),” Soviet Math. Dokl. 27, 372–376 (1983).

22. P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forward-backward splitting,” Multiscale Model Simul. 4, 1168–1200 (2006). [CrossRef]  

23. Y. Nesterov, “Gradient methods for minimizing composite functions,” Math. Prog. 140, 125–161 (2013). [CrossRef]  

24. H. Attouch, J. Bolte, and B. F. Svaiter, “Convergence of descent methods for semi-algebraic and tame problems: Proximal algorithms, forward–backward splitting, and regularized Gauss–Seidel methods,” Math. Prog. 137, 91–129 (2013). [CrossRef]  

25. A. Beck and M. Teboulle, “Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems,” IEEE Trans. Image Process. 18, 2419–2434 (2009). [CrossRef]   [PubMed]  

26. A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” J. Math. Imaging Vision 40, 120–145 (2011). [CrossRef]  

27. U. S. Kamilov, “A parallel proximal algorithm for anisotropic total variation minimization,” IEEE Trans. Image Process. 26, 539–548 (2017). [CrossRef]  

28. S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Found. Trends Mach. Learn. 3, 1–122 (2011). [CrossRef]  

29. S. Lefkimmiatis, J. Ward, and M. Unser, “Hessian Schatten-norm regularization for linear inverse problems,” IEEE Trans. Image Process. 22, 1873–1888 (2013). [CrossRef]   [PubMed]  

30. L. Bottou, “Large-scale machine learning with stochastic gradient descent,” in “Proceedings of COMPSTAT’2010: 19th International Conference on Computational StatisticsParis France, August 22–27, 2010 Keynote, Invited and Contributed Papers,” Y. Lechevallier and G. Saporta, eds. (Physica-Verlag HD, Heidelberg, 2010), pp. 177–186.

31. A. J. Devaney, Mathematical Foundations of Imaging, Tomography and Wavefield Inversion (Cambridge University Press, 2012). [CrossRef]  

32. J. A. Stratton, Electromagnetic Theory (John Wiley & Sons, 2007).

33. M. Unser, E. Soubies, F. Soulez, M. McCann, and L. Donati, “GlobalBioIm: A unifying computational framework for solving inverse problems,” in “Proceedings of the OSA Imaging and Applied Optics Congress on Computational Optical Sensing and Imaging (COSI’17),” (San Francisco CA, USA, 2017). Paper no. CTu1B.

34. Jean-Michel Geffrin, Pierre Sabouroux, and Christelle Eyraud, “Free space experimental scattering database continuation: experimental set-up and measurement precision,” Inverse Probl. 21, 6(2005). [CrossRef]  

35. J.-J. Moreau, “Fonctions convexes duales et points proximaux dans un espace hilbertien,” C. R. Acad. Sci. Ser. A Math. 255, 2897–2899 (1962).