Versatile reconstruction framework for diffraction tomography with intensity measurements and multiple scattering

Thanh-An Pham; Emmanuel Soubies; Alexandre Goy; Joowon Lim; Ferréol Soulez; Demetri Psaltis; Michael Unser

doi:10.1364/OE.26.002749

1. Introduction

Inverse scattering finds a large range of applications going from microwave imaging [1, 2] to microscopy [3]. For example, in biology, this technique referred to as optical diffraction tomography (ODT) is a method of choice to observe transparent samples without staining [4]. Given complex measurements of the scattered field when the sample is illuminated with tilted incident waves, it allows one to do quantitative imaging of the distribution of refractive indices (RI) through the resolution of an inverse scattering problem [5–8] (Section 1.1.1). Because in some applications the phase of the scattered field cannot be measured accurately, alternative methods rely on intensity measurements only (Section 1.1.2). In microscopy, this allows to overcome the interferometric system required to record holographic data. The price to pay, however, is that the reconstruction problem becomes more challenging. In practice, it is usually addressed by alternating between phase retrieval and RI estimation [9–12].

1.1. Related works

1.1.1. Reconstruction from complex measurements

In the early works [13, 14], RI reconstruction is performed using a suitable variant of filtered back-projection algorithm (FBP) [15]. This kind of approach is computationally efficient, but the underlying model ignores the effect of diffraction. Under the assumption that the scattered field is weak compared to the incident one, improved methods rely on the first Born model [4, 16] or the first Rytov model [17, 18]. Both are linearizations of the Lippmann-Schwinger equation governing the ODT acquisition process. The main difference is that the Born model works directly with the scattered field whereas the Rytov model uses the unwrapped phase of the measurements. Several studies have shown that the Rytov model yields more accurate reconstructions than the Born model [18–20]. Yet, the validity of these linear models is restricted to weakly scattering samples.

To overcome this limitation, nonlinear models that account for multiple scattering have been proposed, such as the beam-propagation method (BPM) [7, 21], the contrast source-inversion method (CSI) [1, 2], hybrid methods [6, 22], the conjugate-gradient method (CGM) [23, 24], or the recursive Born approximation [25]. Finally, within a regularized variational approach, iterative forward models that solve the Lippmann-Schwinger equation have been recently used in [5, 8, 22]. This equation allows for a finer description of wave-optics and leads to improved reconstruction results. Note that it is also closely related to the discrete dipole approximation (DDA) or the methods of moments (MoM) used in microwave imaging [26]. A review on ODT reconstruction methods for bioimaging can be found in [3].

1.1.2. RI reconstruction from intensity-only measurements

The more challenging problem of reconstructing refraction indices from intensity-only measurements was recently addressed in the context of Fourier ptychographic microscopy [9, 27–30], although related methods were previously proposed in other domains [31–40]. The conventional reconstruction scheme consists in alternating between recovering the measurement phase and reconstructing the RI of the sample. The phase retrieval is generally performed with the Gerchberg-Saxton (GS) projection operator [41], while the RI reconstruction step is essentially the same as in ODT (intensity and phase). Over the years, the forward models have evolved from linear [9] to nonlinear [10–12]. For more details, we refer the reader to [42]. In Table 1, we have sketched the evolution of reconstruction methods for both holographic (i.e., intensity and phase) and intensity-only measurements.

Table 1. Some RI reconstruction algorithms from holographic (i.e., complex) or intensity-only measurements. Ref.: Reference. Algo.: Algorithm. Reg.: Regularization. rec. Born: Recursive Born. BPM: Beam-propagation method. LSm: Lippmann-Schwinger model. GS: Gerchberg-Saxton projection operator. E.: Embedded. TV: Total-variation constraint. a.h.: ad hoc. 3PIE: Ptychographical iterative engine. GD: Gradient descent. FBS: Forward-backward splitting. iter: Iterative.

View Table | View all tables in this article

1.2. Contributions

In this paper, we leverage recent advances in phase retrieval, nonlinear physical models, and modern regularization. We propose a unified framework that can cope with forward models at various levels of sophistication (e.g., Born [4], beam-propagation method [44], Lippmann-Schwinger [45]) and with various sparse regularizers (e.g., TV [46], Hessian-based [47]). This is possible because of the modularity of the proposed approach, which comes from an adequate splitting of the initial problem into simpler subproblems. Moreover, our method can be easily adapted to different types of noise by the way of specific data-fidelity terms for which an explicit expression of the proximity operator is available. Finally, we validate the proposed method on several simulated and real datasets using both the beam-propagation method (BPM) and the Lippmann-Schwinger forward model (LSm) together with a total-variation (TV) regularizer.

1.3. Organization of the paper

In Section 2, we describe our unified regularized reconstruction framework for the intensity-only problem and develop our optimization strategy. In Section 3, we present the nonlinear forward models used in this work. Finally, in Section 4, we present our numerical results on simulated and real data.

2. Unified regularized reconstruction framework

We consider a 2D spatial domain Ω discretized into M = M_x × M_z pixels with steps δ_x and δ_z. Each cell is characterized by its refractive value n_m and gathered into the vector n ∈ ℝ^M. In Ω, a sample of unknown RI distribution is immersed in a background medium of RI $\tilde{n}$ . We introduce the scattering potential f ∈ ℝ^M, the components of which being defined as $f_{m} = {\tilde{k}}^{2} (n_{m}^{2} / {\tilde{n}}^{2} - 1) \forall m \in [1 \dots M]$ with $\tilde{k} = \frac{2 π \tilde{n}}{λ}$ the background wavenumber. This sample is illuminated by a series of P ∈ ℕ tilted plane waves ${u_{p}^{in} \in ℂ^{M}}_{p \in [1 \dots P]}$ of free-space wavelength λ. The interaction of an incident wave with the sample produces a scattered wave. The sum of the incident and scattered waves yields the total field u_p ∈ ℂ^M. It is worth noting that an alternative setting considers a fixed source-detector relationship while rotating the sample. The general framework that we describe in this work is applicable to both setups.

We represent the ODT forward model by the operator $S : ℝ^{M} \times ℂ^{M} \to ℂ^{N}$ ℝ^M × ℂ^M → ℂ^N. Given the scattering potential f ∈ ℝ^M and an incident wave uⁱⁿ, $S (f, u^{in})$ gives the total field on the detector region Γ (see Fig. 1). The intensity-only measurements {y_p ∈ ℝ^N}_p_∈[1…_P_] are related to the model by

y_{p} = {| S (f, u_{p}^{in}) |}^{2} + η_{p}, \forall p \in [1 \dots P],

where η_p ∈ ℝ^N, ∀p ∈ [1 … P], represent noise and | · |² is a component-wise squared magnitude. Since

u_{p}^{in}

is fixed and known,

S (f, u_{p}^{in})

will be denoted by

S_{p} (f)

thereafter. Note that

S

can be any model of wave scattering. Two specific models will be detailed in Section 3 and used in the experiments of Section 4.

Fig. 1 Optical diffraction tomography setup (intensity-only). A sample with the refractive index n ∈ ℝ^M is immersed in a background medium of index ñ and impinged by an incident plane wave with a given orientation (wave vector k_b). The interaction of the wave with the object produces a scattered wave (forward and backward). The squared magnitude of the total field, which corresponds to the sum of the the incident and scattered waves, is recorded by the detector.

Download Full Size | PDF

Formulation of the inverse problem

Within the context of variational approaches for inverse problems, it is customary to estimate the scattering potential f ∈ ℝ^M from the measurements {y_p ∈ ℝ^N}_p_∈[1…_P_] by solving the optimization problem

\hat{f} \in {\underset{f \in B}{\arg \min} (\sum_{p = 1}^{P} D ({| S_{p} (f) |}^{2}, y_{p}) + τ ℛ (Lf))} .

The functional

D : ℝ^{N} \times ℝ^{N} \to ℝ_{\geq 0}

measures the fidelity of the model to the data. The regularization function

ℛ : ℝ^{K} \to ℝ_{\geq 0}

promotes the sparsity of the quantity Lf, where L: ℝ^M → ℝ^K is a linear operator (e.g., gradient, Hessian). The scalar τ > 0 is a tradeoff parameter that balances the effect of these two terms. The set B represents physical constraints on the scattering potential (e.g., nonnegativity constraint

B = ℝ_{\geq 0}^{M}

). From a Bayesian point of view, we can relate

D

to the log-likelihood of the noise distribution. Because the number of measurements N is much smaller than the number of unknowns M, the data-fidelity term

D

does not generally admit a unique global minimizer. The regularization term ℛ(L·) and the set B should thus be chosen in order to discriminate between candidate solutions using the knowledge that one has on the observed sample.

2.1. Splitting strategy

Inspired by the success of the alternating-direction method of multipliers (ADMM) [48], we propose to split the optimization task in a way that decouples the complex-field-based reconstruction from the phase retrieval. To that end, we introduce the auxiliary variables v_p ∊ ℂ^N, p ∊ [1 … P], and reformulate the problem in Eq. (2) as

(\hat{f}, \hat{v_{1}}, \dots, \hat{v_{P}}) \in {\arg \min_{f, v_{1}, \dots, v_{P}) \in X} (\sum_{p = 1}^{P} D ({| v_{p} |}^{2}, y_{p}) + τ ℛ (Lf))},

where

X = {(f, v_{1}, \dots, v_{P}) \in B \times ℂ^{N \times P} s . t . v_{p} = X_{p} (f) \forall p \in [1 \dots P]} .

The augmented-Lagrangian form of this problem is

ℒ (f, v_{1}, \dots, v_{P}, w_{1}, \dots, w_{P},) = \sum_{p = 1}^{P} D ({| v_{p} |}^{2}, y_{p}) + \frac{ρ}{2} {‖ S_{p} (f) - v_{p} + w_{p} / ρ ‖}_{2}^{2} + τ ℛ (Lf),

where w_p are the Lagrangian multipliers and ρ is a positive scalar. Then, Eq. (5) is minimized using ADMM, which results in the procedure given in Algorithm 1. The problem is now reduced to three simpler subproblems: a phase retrieval that requires the computation of the proximity operator of

D ({| \cdot |}^{2}, y_{p})

, an RI reconstruction problem from complex measurements, and the Lagrangian update of w_p.

Algorithm 1. ADMM for solving Eq. (5)

View Table | View all tables in this article

2.2. Proximity operator for phase-retrieval

At Step 5 of Algorithm 1, we must compute the proximity operator of $\frac{1}{ρ} D ({| \cdot |}^{2}, y_{p})$ , like in

{prox}_{\frac{1}{ρ} D ({| \cdot |}^{2}, y_{p})} (x) = \arg \min_{v \in ℂ^{N}} (\frac{1}{2} {‖ v - x ‖}_{2}^{2} \frac{1}{ρ} D ({| v |}^{2}, y_{p})) .

Here, we benefit from the closed-form expressions that have been recently derived for Gaussian likelihood in [49]. In the present work, we consider the weighted quadratic data-fidelity term

D ({| v |}^{2}, y_{p}) = \frac{1}{2} {‖ {| v |}^{2} - y_{p} ‖}_{w_{p}}^{2},

where

W_{p} = d i a g (w_{1}^{p}, \dots, w_{N}^{p})

is a diagonal matrix and

{‖ \cdot ‖}_{W_{p}}

a weighted ℓ₂-norm such that

{‖ v ‖}_{W_{p}}^{2} = v^{T} W_{p} v

. This scheme can be tuned to two scenarios.

Log-likelihood for Gaussian measurement noise: We set $w_{n}^{p} = 1 / σ^{2}$ , where σ² is the variance of the noise.
Log-likelihood for Poisson measurement noise: We set $w_{n}^{p} = \max {(y_{p, n}, b)}^{- 1}$ , where the minimal value b > 0 accounts for background emission and the dark current of the detector

Following [49], the proximity operator of

D ({| \cdot |}^{2}, y_{p})

given by Eq. (7) is computed component-wise according to

\forall x \in ℂ^{N}, {[{prox}_{\frac{1}{ρ} D ({| \cdot |}^{2}, y_{p})} (x)]}_{n} = ϱ_{n} e^{i \arg (x_{n})},

where ϱ_n is the positive root of the cubic polynomial in ϱ

q_{G} (ϱ) = \frac{4 w_{n}^{p}}{ρ} ϱ^{3} + ϱ (1 - \frac{4 w_{n}^{p}}{ρ} y_{p, n}) - | x_{n} |

which can be efficiently found with Cardano’s method. Note that a closed-form expression has also been derived for an exact model of Poisson noise [49].

2.3. Reconstruction from complex fields

The minimization over f (Step 6 of Algorithm 1) can be achieved by deploying an accelerated forward-backward splitting (FBS) algorithm [50–52]. Two quantities needs to be computed.

The proximity operator of the regularization term $f \mapsto τ / ρ ℛ (Lf)$ , which reads [53] τ ${prox}_{(\frac{τ}{ρ}) ℛ (L \cdot)} (f) = \arg \min_{g \in ℝ_{\geq 0}^{M}} (\frac{1}{2} {‖ g - f ‖}_{2}^{2} + (\frac{τ}{ρ}) ℛ (Lg)) .$ Depending on the choice of $ℛ$ and L, Eq. (10) can admit a closed-form solution or be computed using efficient convex-optimization algorithms (e.g., ADMM [48], Chambolle-Pock [54]).
The gradient of the quadratic term $ℋ : f \mapsto \frac{ρ}{2} {\sum_{p = 1}^{P} ‖ S_{p} (f) - z_{p}^{(k + 1)} ‖}_{2}^{2}$ with $z_{p}^{(k + 1)} = (v_{p}^{(k + 1)} - w_{p}^{(k)} / ρ)$ which, using classical differential rules, is given by $\nabla ℛ (f) = ρ \sum_{p = 1}^{P} Re (J_{S_{p}}^{H} (f) (S_{p} (f) - z_{p}^{(k + 1)})),$ where $J_{S_{p}} (f)$ denotes the Jacobian matrix of $S_{p}$ at f. Hence, given a forward model $S_{p}$ , one needs to provide an efficient computation of the Jacobian matrix $J_{S_{p}} (f)$ .

The classical FBS algorithm has been shown to converge in the nonconvex case, albeit locally only [55]. To the best of our knowledge, there exists no theoretical proof of convergence of its accelerated version for nonconvex cases. However, we observed that the accelerated FBS always converged in our experiments.

3. Solutions for specific nonlinear forward models

While the reconstruction algorithm developed in Section 2 can deal with any forward model ( $S (f)$ in Eq. (1)), we shall focus here on two specific nonlinear models that are known to outperform linear models [5, 7, 8]. For the sake of clarity, we drop the illumination index p and present the forward model for a generic incident wave uⁱⁿ ∈ ℂ^M.

3.1. Beam-propagation method

BPM computes the total field slice-by-slice along the optical axis z [44, 56]. In the following, indexing a vector by z refers to its restriction to the corresponding z-slice. Given an incident wave uⁱⁿ, the total field u on Ω is computed recursively according to

{\begin{array}{l} u_{0} (f) = u_{0}^{in} \\ u_{z^{'}} (f) = (u_{z^{'} - δ_{z}} (f) * h_{prop}^{δ_{z}}) ⊙ p_{z^{'}} (f) \end{array}

where denotes the Hadamard product, the convolution operation, and δ_z the step size in the axial direction. In Eq. (12),

u_{z^{'} - δ_{z}} (f)

is first propagated to the axial position z′ by convolution with the propagation kernel

h_{prop}^{δ_{z}} = ℱ^{- 1} {\exp (i δ_{z} (\sqrt{{\tilde{k}}^{2} 1_{M_{x}} - k_{x}^{2}}))} (diffraction step),

where ℱ⁻¹ is the 1D discrete inverse Fourier transform and

k_{x} \in ℝ^{M_{x}}

is the frequency variable along the x direction. This convolution is followed by a point-wise multiplication with the phase mask

p_{z^{'}} (f) = \exp (i {\tilde{k}}_{0} δ_{z} δ n_{z^{'}} (f)) (refraction step) .

Here, we introduce the RI variation δn(f) ∈ ℝ^M which is related to f through

δ n (f) = \tilde{n} (\sqrt{1 + \frac{f}{\tilde{k}} - 1}) .

Once the axial position z′ reaches the last slice L_z = M_zδ_z of Ω,

\tilde{G} : ℂ^{M_{x}} \to ℂ^{N}

propagates the total field to the detector plane and restricts it to the detector positions. The BPM forward model is defined by

S (f) = \tilde{G} (u_{L_{z}} (f)) .

Note that the BPM forward model can also be directly defined with δn [7].

Computation of the Jacobian matrix

The recursive structure of the BPM forward model enables an efficient computation of the Jacobian $J_{S} (f)$ using an error-backpropagation scheme [7].

3.2. Iterative Lippmann-Schwinger model

In its discrete form, the Lippmann-Schwinger equation reads

u = u^{in} + G d i a g (f) u,

where G ∈ ℂ^M^×^M denotes the matrix defining the convolution (over Ω) with the Green’s function

G (r) = {\begin{array}{l} \frac{i}{4} H_{0}^{(1)} (\tilde{k} {‖ r ‖}_{2}), & in 2 D \\ \frac{1}{4 π} \frac{e^{i \tilde{k} {‖ r ‖}_{2}}}{{‖ r ‖}_{2}} & in 3 D . \end{array}

From Eq. (17), the total field on Ω is computed by solving the quadratic-minimization problem

u (f) = \arg \min_{u \in ℂ^{M}} (\frac{1}{2} {‖ (I - G d i a g (f)) u - u^{in} ‖}_{2}^{2})

using an efficient algorithm such as the conjugate gradient (CG). Based on Eq. (19), the Lippmann-Schwinger iterative forward model (LSm) is defined by

S (f) = {\tilde{G} d i a g (f) u (f) + u^{in} |}_{Γ},

where

\tilde{G}

encodes the convolution with the Green’s function evaluated at the sensor positions Γ and uⁱⁿ|_Γ denotes the incident field uⁱⁿ on the area Γ.

Computation of the Jacobian matrix

As shown in [5], the associated Jacobian matrix $J_{S} (f)$ is given by

J_{S} (f) = \tilde{G} (I + d i a g (f) {(I - G d i a g (f))}^{- 1} G) d i a g (u (f))

and can also be estimated using the CG algorithm, if desired.

4. Results

We first assess the suitability of our framework to reconstruct simulated samples. Then, we validate our approach on experimental data. Finally, we evaluate the performance of the method for limited measurements.

We compare the solutions of our framework to those obtained with the light field refocusing (LFR) method [57] which were also used as initial guesses for Algorithm 1. For the regularizer $ℛ (L \cdot)$ , we use the TV norm, known to attenuate the missing-cone problem. Moreover, we enforce a nonnegativity constraint on the scattering potential by setting $B = ℝ_{\geq 0}^{M}$ . Because the RI reconstruction step can be computationally intensive, we adopted acceleration strategies. The gradient in Eq. (11) was computed for a subset of the angles [1 … P]. This subset was changed at each iteration while keeping a constant angle difference between them. We implemented the algorithms using an inverse-problem library developed in our group [58] (GlobalBioIm: http://bigwww.epfl.ch/algorithms/globalbioim/).

4.1. Validation on simulated data

Simulation setup

We consider the three samples presented in Fig. 2 (top row). They are immersed in water $(\tilde{n} = 1.33)$ as well as the source and the sensors. They were impinged by thirty-one incident waves with angles ranging from −45° to 45°. These waves were propagated from the bottom to the top of the (33λ × 33λ domain with λ = 406 nm. Simulations were performed on a fine grid (1024 × 1024) with a pixel area of (0.03λ)² using the LSm forward model. The 1024 sensors are evenly placed on a straight line of length 33λ above the sample at 16.5λ from the center. Finally, these measurements were reduced to N = 512 using averaging.

Fig. 2 RI distribution for three fibers, a simulated cell, and the Shepp-Logan in the first, second, and third column, respectively. The ground truth and the reconstructions from the LFR, BPM, and LSm proposed methods are shown in Row 1 to 4, respectively. The samples are immersed in water $(\tilde{n} = 1.33)$ . Thirty-one views were acquired with a tilted plane-wave illumination. The angles ranged from −45° to 45°. The sample is illuminated from below. The 1024 sensors are evenly placed on a straight line of length 33λ above the sample at 16.5λ from the center. The measurements were reduced to 512 using averaging.

Download Full Size | PDF

Reconstruction parameters

The regularization parameter τ was selected in order to minimize the relative error ${‖ \hat{x} - x_{true} ‖}_{2} / {‖ x_{true} ‖}_{2}$ . The outer loop in Algorithm 1 (inner FBS at step 6 in Algorithm 1, respectively) was stopped when the relative change between two iterates is below 10⁻⁸ or after 20 (50, respectively) iterations. The step size in the FBS algorithm was set to 5 · 10⁻⁴ and 5 · 10⁻³ for BPM and LSm, respectively. The penalty parameter was set to ρ = 10⁻³. Finally, the reconstructions were performed on a (512 × 512) grid. The regularization parameter τ was tuned by hand for each sample.

Discussion

As shown in Fig. 2, the proposed framework is able to reconstruct the samples despite the lack of phase information. Both forward models obviously perform better than LFR which only relies on geometrical optics. We observe that the LSm forward model yields better reconstructions than the BPM forward model. Our framework with LSm is able to retrieve most details of the object. The shape as well as the RI of the samples are well estimated. These observations are quantified by the relative error presented in Table 2.

Table 2. Reconstruction performance. The relative error $ϵ = \frac{{‖ \hat{x} - x_{true} ‖}_{2}}{{‖ x_{true} ‖}_{2}}$ is shown. The proposed method with BPM was 3 to 6 times faster than with LSm.

View Table | View all tables in this article

4.2. Validation on experimental data

We validate our framework using the publicly available experimental datasets of the Institut Fresnel [59]. In this experiment, the sensors were circularly distributed around the objects at a distance of 1.67 m from their centers with 1° step (see Fig. 3). Eight sources (E₁₋₈), uniformly distributed around the object, were activated sequentially. At each activation, the total field was measured using 241 sensors (S₂₄₁), excluding the sensors closer than 60° from the source. We reconstructed the three targets FoamDielExt, FoamDielInt, and FoamTwinDiel using the TM polarization at 3 GHz (i.e., λ = 10 cm). Each two-dimensional inhomogeneous sample is depicted in Fig. 4 (top row). The indicated permittivities were experimentally measured and are subject to uncertainties [59].

Fig. 3 Acquisition setup for the Fresnel dataset. The sensors (dots in the inner circle) correspond to the illumination angle of 0° (i.e., E₁). The measurements are restricted either by reducing the number of sensors (sensors sets S₂₄₁, …, S₉₁) or the number of acquired views (emitters E₁, …, E₈).

Download Full Size | PDF

Fig. 4 Permittivity reconstruction of the Fresnel datasets by LSm. From left to right: FoamDielExt, FoamDielInt, and FoamTwinDiel. From top to bottom: ground truth, reconstructions from complex (using [5]) and intensity-only (proposed method) measurements, respectively, and magnitude and phase of the predicted (solid curve) vs true (dashed curve) measurements (0° illumination angle). The two curves often overlap. For the solutions from complex measurements, the regularization parameters were set at 1.6 · 10⁻², 3 · 10⁻³, and 9 · 10⁻³ for FoamDielExt, FoamDielInt, and FoamTwinDiel, respectively. For the solutions from intensity-only measurements, the regularization parameters were set at 7 · 10⁻², 9 · 10⁻³, and 4 · 10⁻² for FoamDielExt, FoamDielInt, and FoamTwinDiel, respectively.

Download Full Size | PDF

For the reconstruction, we consider a (15 × 15 cm²) area discretized over (256 × 256) pixels. This yields a pixel area of about (0.0586 cm)². We reconstruct these samples with the LSm forward model and compare the results with the RI reconstruction from holographic measurements [5].

The chosen regularization parameters τ were tuned by hand for both algorithms. The penalty parameter ρ was set to 5. We initialized both algorithms with the background value.

The results presented in Fig. 4 suggest that the two methods perform similarly in terms of quality. The shape and the permittivity of the samples are both remarkably well recovered despite the high contrasts. Furthermore, the bottom graphs in Fig. 4 show that the retrieved phase corroborates the measurement data for each sample. The similar performances observed for these samples suggest the intensity measurements still contain some phase information due to the diffraction.

4.3. Reducing the number of measurements

In this section, we assess the effect of a reduction in the number of measurements. To that end, we combined two methods. On one hand, we incrementally ignored illumination angles. On the other hand, we reduced the number of sensors, starting from no restriction (i.e., S₂₄₁) to the smallest set of sensors S₉₁ (see Fig. 3). This strategy allowed us to explore the missing-cone problem. By progressively limiting the available measurements, we converged to a setup similar to that of tomographic phase microscopy [3]. The reconstruction obtained for the easiest scenario (i.e., 8 views and S₂₄₁) was considered as a reference. Then, the regularization parameter τ was tuned in order to minimize the relative error with respect to this reference.

As shown in Fig. 5, the quality of the reconstructions is remarkable, even in extreme cases. This is due to the use of modern regularization.

Fig. 5 Permittivity reconstructions of the Fresnel dataset with a limited number of measurements. From left to right: P = 3, 5, 7, and 8 views were used to reconstruct the sample FoamDielExt. From top to bottom: The sensors were included in the sets S₂₄₁, S₁₈₁, S₁₅₁, S₁₂₁, and S₉₁, respectively. The reconstruction error with respect to the best solution (i.e., E₈, S₂₄₁) is shown at the top left of each image.

Download Full Size | PDF

5. Conclusion

We have proposed a variational formulation of the reconstruction of refractive indices from intensity-only measurements. It allows us to take advantage of efficient algorithms to solve subproblems. Our framework is able to handle several forward models and any regularization. Notably, we showed that the iterative Lippmann-Schwinger model combined with total-variation regularization reconstructs highly scattering samples from intensity-only measurements, even in ill-posed configurations. Furthermore, our results suggest our method can reconstruct RI samples in even more difficult cases where few measurements are available.

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. A. Abubakar and P. van den Berg, “The contrast source inversion method for location and shape reconstructions,” Inverse Probl. 18, 495 (2002). [CrossRef]

2. P. M. van den Berg, A. Van Broekhoven, and A. Abubakar, “Extended contrast source inversion,” Inverse Probl. 15, 1325 (1999). [CrossRef]

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

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

5. E. Soubies, T.-A. Pham, and M. Unser, “Efficient inversion of multiple-scattering model for optical diffraction tomography,” Opt. Express 25, 21786–21800 (2017). [CrossRef] [PubMed]

6. E. Mudry, P. 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]

7. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Optical tomographic image reconstruction based on beam propagation and sparse regularization,” IEEE Trans. Comput. Imaging 2, 59–70 (2016). [CrossRef]

8. H. Y. Liu, D. Liu, H. Mansour, P. T. Boufounos, L. Waller, and U. S. Kamilov, “SEAGLE: Sparsity-driven image reconstruction under multiple scattering,” IEEE Trans. Comput. Imaging, https://arxiv.org/abs/1705.04281 (2017).

9. G. Zheng, R. Horstmeyer, and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nat. Photonics 7, 739–745 (2013). [CrossRef]

10. A. M. Maiden, M. J. Humphry, and J. M. Rodenburg, “Ptychographic transmission microscopy in three dimensions using a multi-slice approach,” J. Opt. Soc. Am. A 29, 1606–1614 (2012). [CrossRef]

11. L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an LED array microscope,” Optica 2, 104–111 (2015). [CrossRef]

12. P. Li, D. J. Batey, T. B. Edo, and J. M. Rodenburg, “Separation of three-dimensional scattering effects in tilt-series Fourier ptychography,” Ultramicroscopy 158, 1–7 (2015). [CrossRef] [PubMed]

13. F. Charrière, A. Marian, F. Montfort, J. Kuehn, T. Colomb, E. Cuche, P. Marquet, and C. Depeursinge, “Cell refractive index tomography by digital holographic microscopy,” Opt. Lett. 31, 178–180 (2006). [CrossRef] [PubMed]

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

15. A. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (SIAM, 2001), Vol. 33. [CrossRef]

16. M. Debailleul, B. Simon, V. Georges, O. Haeberlé, and V. Lauer, “Holographic microscopy and diffractive microtomography of transparent samples,” Meas. Sci. Technol. 19, 074009 (2008). [CrossRef]

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

18. 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]

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

20. J. Lim, K. Lee, K. 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]

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

22. T. Zhang, C. Godavarthi, P. C. Chaumet, G. Maire, H. Giovannini, A. Talneau, M. Allain, K. Belkebir, and A. Sentenac, “Far-field diffraction microscopy at λ/10 resolution,” Optica 3, 609–612 (2016). [CrossRef]

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

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

25. 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]

26. M. Yurkin and A. Hoekstra, “The discrete dipole approximation: An overview and recent developments,” J. Quant. Spectrosc. Radiat. Transf. 106, 558–589 (2007). [CrossRef]

27. L. Tian, X. Li, K. Ramchandran, and L. Waller, “Multiplexed coded illumination for Fourier ptychography with an LED array microscope,” Biomed. Opt. Express 5, 2376–2389 (2014). [CrossRef] [PubMed]

28. L.-H. Yeh, J. Dong, J. Zhong, L. Tian, M. Chen, G. Tang, M. Soltanolkotabi, and L. Waller, “Experimental robustness of Fourier ptychography phase retrieval algorithms,” Opt. Express 23, 33214–33240 (2015). [CrossRef]

29. Y. Zhang, W. Jiang, L. Tian, L. Waller, and Q. Dai, “Self-learning based Fourier ptychographic microscopy,” Opt. Express 23, 18471–18486 (2015). [CrossRef] [PubMed]

30. L. Bian, J. Suo, G. Zheng, K. Guo, F. Chen, and Q. Dai, “Fourier ptychographic reconstruction using Wirtinger flow optimization,” Opt. Express 23, 4856–4866 (2015). [CrossRef] [PubMed]

31. G. Gbur and E. Wolf, “The information content of the scattered intensity in diffraction tomography,” Inf. Sci. 162, 3–20 (2004). [CrossRef]

32. T. C. Wedberg and J. J. Stamnes, “Experimental examination of the quantitative imaging properties of optical diffraction tomography,” J. Opt. Soc. Am. A 12, 493–500 (1995). [CrossRef]

33. M. Holler, M. Guizar-Sicairos, E. H. R. Tsai, R. Dinapoli, E. Muller, O. Bunk, J. Raabe, and G. Aeppli, “High-resolution non-destructive three-dimensional imaging of integrated circuits,” Nature 543, 402–406 (2017). [CrossRef] [PubMed]

34. M. H. Maleki, A. J. Devaney, and A. Schatzberg, “Tomographic reconstruction from optical scattered intensities,” J. Opt. Soc. Am. A 9, 1356–1363 (1992). [CrossRef]

35. M. H. Maleki and A. J. Devaney, “Phase-retrieval and intensity-only reconstruction algorithms for optical diffraction tomography,” J. Opt. Soc. Am. A 10, 1086–1092 (1993). [CrossRef]

36. T. Takenaka, D. J. Wall, H. Harada, and M. Tanaka, “Reconstruction algorithm of the refractive index of a cylindrical object from the intensity measurements of the total field,” Microw. Opt. Technol. Lett. 14, 182–188 (1997). [CrossRef]

37. M. Lambert and D. Lesselier, “Binary-constrained inversion of a buried cylindrical obstacle from complete and phaseless magnetic fields,” Inverse Probl. 16, 563 (2000). [CrossRef]

38. L. Crocco, M. D’Urso, and T. Isernia, “Inverse scattering from phaseless measurements of the total field on a closed curve,” J. Opt. Soc. Am. A 21, 622–631 (2004). [CrossRef]

39. A. Litman and K. Belkebir, “Two-dimensional inverse profiling problem using phaseless data,” J. Opt. Soc. Am. A 23, 2737–2746 (2006). [CrossRef]

40. W. Van den Broek and C. T. Koch, “Method for retrieval of the three-dimensional object potential by inversion of dynamical electron scattering,” Phys. Rev. Lett. 109, 245502 (2012). [CrossRef]

41. R. Gerchberg and W. O. Saxton, “Phase determination from image and diffraction plane pictures in electron-microscope,” Optik 34, 275–284 (1971).

42. K. Guo, S. Dong, and G. Zheng, “Fourier ptychography for brightfield, phase, darkfield, reflective, multi-slice, and fluorescence imaging,” IEEE J. Sel. Topics Quantum Electron. 22, 77–88 (2016). [CrossRef]

43. D. Ren, E. Bostan, L. Yeh, and L. Waller, “Total-variation regularized Fourier ptychographic microscopy with multiplexed coded illumination,” in “Mathematics in Imaging,” (Optical Society of America, 2017), paper MM3C–5. [CrossRef]

44. J. Van Roey, J. Van der Donk, and P. Lagasse, “Beam-propagation method: Analysis and assessment,” J. Opt. Soc. Am. A 71, 803–810 (1981). [CrossRef]

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

46. L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D 60, 259–268 (1992). [CrossRef]

47. 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]

48. 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]

49. F. Soulez, E. Thiébaut, A. Schutz, A. Ferrari, F. Courbin, and M. Unser, “Proximity operators for phase retrieval,” Appl. Opt. 55, 7412–7421 (2016). [CrossRef] [PubMed]

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

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

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

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

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

55. H. Attouch, J. Bolte, and B. 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]

56. J. Fleck, J. Morris, and M. Feit, “Time-dependent propagation of high energy laser beams through the atmosphere,” Appl. Phys. A 10, 129–160 (1976). [CrossRef]

57. G. Zheng, C. Kolner, and C. Yang, “Microscopy refocusing and dark-field imaging by using a simple LED array,” Opt. Lett. 36, 3987–3989 (2011). [CrossRef] [PubMed]

58. 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 CTu1B.

59. J. Geffrin, P. Sabouroux, and C. Eyraud, “Free space experimental scattering database continuation: experimental set-up and measurement precision,” Inverse Probl. 21, S117 (2005). [CrossRef]

Require: {y_p}_p_∈[1…_P_], $f \in ℝ_{\geq 0}^{M}$ , ρ > 0, τ > 0
1:	f⁽⁰⁾ = f
2:	$w_{p}^{(0)} = 0_{ℂ^{N}}$ , ∀p ∈ [1 … P]
3:	k = 0
4:	while (not converged) do
5:	$v_{p}^{(k + 1)} = {prox}_{\frac{1}{ρ} D ({\| \cdot \|}^{2}, y_{p})} (S_{p} (f^{(k)}) + \frac{w_{p}^{(k)}}{ρ})$ , ∀p ∈ [1 … P] {Phase retrieval}
6:	$f^{(k + 1)} = \arg \min_{f \in B} (\frac{ρ}{2} \sum_{p = 1}^{P} {‖ S_{p} (f) - v_{p}^{(k + 1)} + w_{p}^{(k)} / ρ ‖}_{2}^{2} + τ ℛ (Lf))$ {RI reconstruction}
7:	$w_{p}^{(k + 1)} = w_{p}^{(k)} + ρ (S_{p} (f^{(k + 1)}) - v_{p}^{(k + 1)})$ ∀p ∈ [1 … P] {Update Lagrangian}
8:	k ← k + 1
9:	end while
10:	return f⁽^k⁾

ϵ	Three fibers	Cell	Shepp-Logan
LFR	1.4 · 10⁻²	1.36 · 10⁻²	2.06 · 10⁻²
BPM	4.74 · 10⁻³	6.05 · 10⁻³	1.28 · 10⁻²
LSm	1.33 · 10⁻³	4.04 · 10⁻³	1.02 · 10⁻²

Require: {y_p}_p_∈[1…_P_], $f \in ℝ_{\geq 0}^{M}$ , ρ > 0, τ > 0
1:	f⁽⁰⁾ = f
2:	$w_{p}^{(0)} = 0_{ℂ^{N}}$ , ∀p ∈ [1 … P]
3:	k = 0
4:	while (not converged) do
5:	$v_{p}^{(k + 1)} = {prox}_{\frac{1}{ρ} D ({\| \cdot \|}^{2}, y_{p})} (S_{p} (f^{(k)}) + \frac{w_{p}^{(k)}}{ρ})$ , ∀p ∈ [1 … P] {Phase retrieval}
6:	$f^{(k + 1)} = \arg \min_{f \in B} (\frac{ρ}{2} \sum_{p = 1}^{P} {‖ S_{p} (f) - v_{p}^{(k + 1)} + w_{p}^{(k)} / ρ ‖}_{2}^{2} + τ ℛ (Lf))$ {RI reconstruction}
7:	$w_{p}^{(k + 1)} = w_{p}^{(k)} + ρ (S_{p} (f^{(k + 1)}) - v_{p}^{(k + 1)})$ ∀p ∈ [1 … P] {Update Lagrangian}
8:	k ← k + 1
9:	end while
10:	return f⁽^k⁾

ϵ	Three fibers	Cell	Shepp-Logan
LFR	1.4 · 10⁻²	1.36 · 10⁻²	2.06 · 10⁻²
BPM	4.74 · 10⁻³	6.05 · 10⁻³	1.28 · 10⁻²
LSm	1.33 · 10⁻³	4.04 · 10⁻³	1.02 · 10⁻²

Versatile reconstruction framework for diffraction tomography with intensity measurements and multiple scattering

Abstract

1. Introduction

1.1. Related works

1.1.1. Reconstruction from complex measurements

1.1.2. RI reconstruction from intensity-only measurements

1.2. Contributions

1.3. Organization of the paper

2. Unified regularized reconstruction framework

Formulation of the inverse problem

2.1. Splitting strategy

2.2. Proximity operator for phase-retrieval

2.3. Reconstruction from complex fields

3. Solutions for specific nonlinear forward models

3.1. Beam-propagation method

Computation of the Jacobian matrix

3.2. Iterative Lippmann-Schwinger model

Computation of the Jacobian matrix

4. Results

4.1. Validation on simulated data

Simulation setup

Reconstruction parameters

Discussion

4.2. Validation on experimental data

4.3. Reducing the number of measurements

5. Conclusion

Funding

Acknowledgments

References and links

Cited By

Figures (5)

Tables (3)

Equations (21)

Optics Express

Holographic measurements				Intensity-only measurements

Ref.	Model	Algo.	Reg.	Ref.	Model	Phase	RI	Reg.
[13,14]	FBP	direct	–						linear
[4,16]	Born	direct	–	[9,27–31]	Born	GS	direct	–
[17,18]	Rytov	direct	–	[34,35]	Rytov	E	iter	–
				[43]	Born	GS	iter	TV
[20]	Rytov	iter	TV
				[10,12]	BPM	GS	3PIE	–
				[11,40]	BPM	GS / E	GD	a.h.
[23,24]	LSm	iter	–	[36,38]	LSm	E / a.h.	iter	–
[22]	LSm	iter	a.h.	[37,39]	LSm	E	iter	a.h.
[7,21]	BPM	FBS	TV	Proposed framework					nonlinear modern regul.
[25]	rec. Born	FBS	TV ↔
[5,8]	LSm	FBS	TV