Comparative study of iterative reconstruction algorithms for missing cone problems in optical diffraction tomography

JooWon Lim; KyeoReh Lee; Kyong Hwan Jin; Seungwoo Shin; SeoEun Lee; YongKeun Park; Jong Chul Ye

doi:10.1364/OE.23.016933

1. Introduction

Measurements of refractive index (RI) provide a means to non-invasively and quantitatively investigate the morphological and chemical information of a sample [1,2], including structures of cells and tissues [3–6], absorptive structures [7–9], biomechanical properties of cell membrane cortex [10–12], dry mass of cells [13, 14], and biochemical information about specific molecules [15–17]. Due to its non-invasiveness and label-free imaging capability, 3-D RI tomography technique has rapidly grown and been actively utilized for the study of biophysics, cell biology, hematology, microbiology, and neuroscience.

In order to experimentally measure 3-D RI distributions of microscopic objects, several quantitative phase imaging (QPI) or digital holographic microscopy techniques have been developed [18–40]. In general, multiple 2-D optical field images of a sample are obtained with the modulation of illumination angle or sample orientation, from which the reconstruction of 3-D RI maps of a sample is performed via tomography algorithms. Experimental schemes include the rotation of a sample [18, 19], the tilting of illumination beams [20–23], the rotation method which rotates the coaxially lined illumination laser source and detector [41], the scanning of illumination wavelengths [24–26], and the axial scanning of a sample illuminated with polychromatic light [27–29].

From multiple 2-D holograms obtained with various illumination angles, a 3-D tomogram is reconstructed either by the projection algorithm or the diffraction algorithm. The diffraction algorithm considers the effect of light diffraction whereas the projection algorithm assumes that the illumination beam propagates in a straight path integrating RI of a sample [21, 22, 42, 43]. Each measured 2-D hologram with various illumination angles is mapped in Fourier space as a plane surface (projection algorithm) or a 2-D hemispheric surface or Ewald sphere (diffraction algorithm) object function distribution, according to the used reconstruction algorithm. Then, the inverse 3-D Fourier transformation of the objective function provides a 3-D RI tomogram of the sample.

Importantly, regardless of experimental schemes and reconstruction algorithms, there exists missing information in the 3-D Fourier spectrum for both the projection and diffraction algorithms due to the limited angle of acceptance of the imaging system determined by the numerical aperture (NA) of an objective lens. Even with a high-NA objective lens, only fraction of diffracted light from a sample can be utilized for the tomographic reconstruction, causing so-called the missing cone problem. The missing cone problem results in the underestimated RI values in reconstructed tomograms and the elongation of the reconstructed shape of a sample along the optical axis [40, 43].

The missing cone (or data) problem is not unique in RI tomography, but has been observed in other imaging modalities such as X-ray computed tomography (CT), electron microscopy (EM), and magnetic resonance imaging (MRI). Addressing the missing cone problem is one of the major research thrusts in these fields [44, 45]. A general approach to address this missing cone problem is to use an iterative algorithm by incorporating a prior information. Gerchberg-Papoulis (GP) algorithm [46, 47] imposes a simple non-negativity constraint to alleviate the under-estimation of RI values. If an unknown image belongs to a function class that consists of piecewise smooth regions separated by sharp edges, edge preserving (EP) regularization can also be employed so that it preserves edges but smooth out noise [48, 49]. A prior knowledge about narrow edges tends to reduce the elongation of a reconstructed image along the optical axis caused by the missing cone problem. Charbonnier et al. proposed a non-convex EP regularization using the half-quadratic regularization technique [48]. Delaney et al. further extended it to the limited angle tomography [49]. Among various EP regularizations, total variation (TV) regularization is one of the most frequently used regularization methods in medical imaging fields such as the compressed sensing reconstruction in MRI [45]. TV regularization is similar to the EP regularization except that TV regularization uses the l₁ norm of the gradient of an image to convexify the problem.

Even the missing cone (data) problem has been extensively addressed in the fields of X-ray CT, MRI, and EM, this issue has not been fully explored in optical RI tomography. For example, TV was theoretically studied in the work by Candes et al, however, it is for random Fourier samples, and the sampling pattern for our ODT problem does not satisfy the RIP conditions in the work [50] which means that extensive verifications are strongly required prior to utilize the algorithm for ODT. Currently, only a few sets of approaches have adapted a specific regularization method in constructing RI tomograms. Sung et al. proposed a regularized recovery scheme for 3-D optical diffraction tomography (ODT) using the non-negativity constraint and EP regularization [40, 51]. TV regularization has been applied to diffraction tomography [52]; a gradient descent method combined with data consistency step was used for the tomographic reconstruction from limited angle data. However, in spite of the importance and implication of the missing problem in RI tomographic reconstruction, a comparative study of the performance of existing regularized inversion algorithms to address the missing cone problem has not been performed in ODT. Here we perform a systematic comparative study of various reconstruction algorithms in the optical RI tomography using both numerical phantoms and experimental data. We implemented and compared three representative algorithms: (i) globally convergent EP algorithm [49] under the first-order Rytov approximation; (ii) TV regularization combined with non-negativity constraint using split Bregman iteration [53, 54]; and (iii) GP algorithm. We discuss the pros and cons of the regularized methods in measuring 3-D RI distributions of microscopic samples.

2. Theory

The principles and details of ODT are well explained in [55]. Therefore, we briefly illustrate the major ideas without extensive descriptions.

2.1. Optical diffraction tomography

In ODT, the fundamental forward equation is given by the scalar Helmholtz equation:

\nabla^{2} U (r) + k^{2} U (r) = - 4 π f (r) U (r),

where

k = \frac{2 π n_{0}}{λ}

, λ is the incident light wavelength in free space, U(r) is a scalar electric field distribution and

f (r) = \frac{k^{2}}{4 π} (\frac{n {(r)}^{2}}{n_{0}^{2}} - 1)

denotes the scattering potential, n(r) is the RI distribution of a sample and n₀ is the RI of a surrounding homogeneous medium. Suppose that the total electrical distribution is given by U(r) = U₀(r) + U_s(r), where U₀(r) is the electric field from the homogeneous surrounding medium and U_s(r) is the field scattered from the sample. Then, under the Born approximation, the scattered field (

U_{s}^{B} (r)

) can be represented as

U_{s}^{B} (r) = \int_{V} f (r^{'}) U_{0} (r^{'}) G (r - r^{'}) d r^{'},

where

G (r - r^{'}) = \frac{e^{i k | r - r^{'} |}}{| r - r^{'} |}

is the Green’s function of 3-D Helmholtz’ equation. For the case of the Rytov approximation, the perturbation is modeled as a phase variation:

U (r) = U_{0} (r) e^{U_{s}^{R}} (r)

, which results in

U_{s}^{R} (r) = \frac{1}{U_{0} (r)} \int_{V} f (r^{'}) U_{0} (r^{'}) G (r - r^{'}) d r^{'} .

2.2. Problem formulation

In this paper, our analysis is based on the first order Rytov approximation. Specifically, we assume that U₀(r) is a plane wave, $U_{0}^{(m)} (r) = e^{i k_{m} \cdot r}$ , where k_m = ks_m is a real-valued wavenumber determined by the directional unit vector for the m-th illumination $s_{m} = (s_{x}^{(m)}, s_{y}^{(m)}, s_{z}^{(m)})$ . Using the Fourier diffraction theorem derived in [56], the measurement of the ODT from the m-th illumination can be written by

g_{m} = \frac{k s_{z}^{(m)}}{2 π i} \int \int_{- \infty}^{+ \infty} U_{0}^{(m)} (r) \ln \frac{U (r)}{U_{0}^{(m)} (r)} e^{i (k s_{x}^{(m)} x + k s_{y}^{(m)} y)} d x d y,

under the assumption that U(r) and

U_{0}^{(m)} (r)

are measured at z = 0 plane. Then, the corresponding forward operator A_m f corresponds to the 3-D Fourier transform of a scattering potential at (k-k_m):

A_{m} f = \int_{V} f (r^{'}) e^{- i (k - k_{m}) \cdot r^{'}} d r^{'} .

Note that Eq. (5) can be calculated using the 2-D Fourier transform after measuring U(r) and $U_{0}^{(m)} (r)$ . Obtaining all 3-D k-space data defined by Eq. (5) for various illumination angles s_m, we can fill in the 3-D k-space of the scattering potential. The corresponding data fidelity term can then be defined as follows

{‖ A f - g ‖}_{2}^{2} = \sum_{m = 1}^{M} ‖ A_{m} f - g_{m} ‖_{2}^{2},

where

A = {[A_{1}^{T} A_{2}^{T} \dots A_{M}^{T}]}^{T}

and

{[g_{1}^{T} g_{2}^{T} \dots g_{M}^{T}]}^{T}

. By adding the regularization and the data fidelity term, we can define the cost function;

J (f) = {‖ A f - g ‖}_{2}^{2} + α R (f),

where R is the regularization function and α is the regularization parameter.

2.3. Iterative reconstruction methods

In this section, we review two typical iterative reconstruction methods that have been used in this study.

2.3.1. Convergent EP regularization with the non-negativity constraint

The cost function for the EP regularization can be defined as following,

J (f) = {‖ A f - g ‖}_{2}^{2} + α \sum_{n} ϕ (D {(f)}_{n}) + β {‖ P f ‖}_{2}^{2},

where ϕ is the potential function that works as an EP regularization,

D (f) = \sqrt{{(\nabla_{x} f)}^{2} + {(\nabla_{y} f)}^{2} + {(\nabla_{z} f)}^{2}}

and P is the diagonal operator that only preserves the positive elements of f. α and β are regularization parameters corresponding to the EP regularization and non-negativity constrain, respectively. Here, ϕ is determined by the prior knowledge about edges satisfying conditions that for

σ (t) = ϕ (\sqrt{t})

, σ: [0,∞) → [0,∞), and is modeled as continuously twice differentiable and strictly concave, with σ (0) = 0,σ′ (0) = 1 and 0 < σ′ (t) ≤ 1 [49]. We used ϕ (t) = T² ln(1 +t/T) [49] satisfying conditions. T is the parameter of σ (t) and it needs to be selected depending on how sharply we want to select edges. Since Eq. (9) is not convex, it is hard to minimize. Therefore, we use the half quadratic regularization [49]. In particular, by plugging an auxiliary variable b, we can rewrite the equation Eq. (9) as [48, 49]

\begin{array}{l} J (f) = \inf_{b} J * (f, b) \\ = \inf_{b} {‖ A f - g ‖}_{2}^{2} + α \sum_{n} (b_{n}) [{(\nabla_{x} f)}_{n}^{2} + {(\nabla_{y} f)}_{n}^{2} + {(\nabla_{z} f)}_{n}^{2} + ψ (b_{n})] + β {‖ P f ‖}_{2}^{2} . \end{array}

where J*(f) is half-quadratic since it is quadratic in f for given b, and also quadratic in b for given f. Therefore we can find the minima of J(f) by alternatingly minimizing with respect to b and f. Here, b^k⁺¹ works as an edge map which can be calculated in the closed form using b^k⁺¹ = σ′(D(f^k)), f^k⁺¹ can be obtained by finding the solution of the following equation,

(A^{H} A + α \nabla_{x}^{T} B^{k + 1} \nabla_{x} + α \nabla_{y}^{T} B^{k + 1} \nabla_{y} + α \nabla_{z}^{T} B^{k + 1} \nabla_{z} + β P^{k + 1}) f^{k + 1} = A^{H} g .

where B^k⁺¹ = diag(b^k+1) and P^k⁺¹ = diag(p^k+1). In this paper, we used the conjugate gradient method to find the solution of Eq. (11). It was shown in [49] that f^k⁺¹ results in convergent subsequences and all their limits are stationary points of Eq. (9).

2.3.2. TV regularization with non-negativity constraint

Consider the constrained minimization problem:

\min_{f} E (f) such that A f = g,

where E(f) is a regularization term related to TV regularization and non-negativity constraint. In TV approach, E(f) = ∑_n D(f) + I(f ≥ 0) where D(f) is the TV term

D (f) = \sqrt{{(\nabla_{x} f)}^{2} + {(\nabla_{y} f)}^{2} + {(\nabla_{z} f)}^{2}}

and I is the indicator function to impose the non-negativity constraint. Since E(f) is not differentiable, it is difficult to minimize Eq. (12). Instead of minimizing the unconstrained cost function Eq. (12), we can modify the problem using Bregman distance and Bregman iteration [57]:

\begin{array}{l} f^{k + 1} = \min_{f} E (f) - < p^{k}, f - f^{k} > + \frac{μ}{2} {‖ A f - g ‖}_{2}^{2}, \\ p^{k}^{+ 1} = p^{k} - μ A^{H} (A f^{k}^{+ 1} - g), \end{array}

By plugging the expression of E(f), we can get

f^{k + 1} = \min_{f} \sum_{n} D (f) + I (f \geq 0) + \frac{μ}{2} {‖ A f - g^{k} ‖}_{2}^{2},

g^{k}^{+ 1} = g^{k} - g - A f^{k}^{+ 1} .

The main technical difficulty lies in solving Eq. (14), which can be addressed using variable splitting. The details can be found in Appendix A.

3. Method

3.1. Experimental setup

Fused silica microsphere (44054, Sigma-Aldrich Inc., USA), red blood cells (RBCs), and hepatocyte cells (Huh-7 cell line, Apath, Brooklyn, NY, USA) were used as experimental samples. All samples were prepared by following the standard protocols [58]. The samples were diluted in Dulbeccos buffered saline (DPBS, RI = 1.337 at 532 nm [59]) and sandwiched between two coverslips before loaded on a microscope (IX-73, Olympus Inc., Japan) which was modified for ODT.

For measuring the 3-D RI tomograms of samples, we implemented tomographic phase microscopy (TPM) based on a Mach-Zehnder interferometer using a laser source (532 nm, 50 mW, Cobolt, Solna, Sweden)(see [9, 43] for details). We used a high-NA objective lens for the illumination (UPLFLN 60×, NA = 0.9 Olympus Inc., Japan) and for the detection (UPLSAPO 60×, NA = 1.42, oil immersion). Holograms were recorded using a CMOS camera (Neo sCMOS, ANDOR Inc., Northern Ireland, UK). More than 300 holograms with different illumination angles were captured, from which a 3-D RI tomogram was retrieved [60]. The illumination angle was scanned in a spiral pattern using a two-axis galvanometric mirror (GVS012/M, Thorlabs, USA). Because the NA of the condenser lens is strongly related to the solid angle of the missing cone, algorithms may show different performances depending on the NA of the used condenser lens. Therefore, we utilized both low- and high-NA cases by controlling the maximum voltage to the galvanometric mirror, which results in the effective NA of NA_eff = 0.5 and 0.8, whose respective theoretical resolutions are 0.14 μm and 0.121 μm, laterally; and 0.525 μm and 0.46 μm, axially, and experimental resolutions are 0.19 μm and 0.14 μm, laterally; and 2.944 μm and 2.706 μm, axially. Here, the axial resolution depends on the lateral structure of the sample due to the gourd bottle shaped ky-kz cross section. Figure 1(a) shows the scanning trajectories for the low (red) and high (green) NA_eff cases, respectively; and Figs. 1(b)–1(c) show the corresponding 3-D k-space sampling trajectories for the illumination beam. To compare numerical and experimental results, microsphere phantoms with the same RI and size of fused silica microsphere were created, and simulated using the experimentally measured angular scanning trajectories. Biological samples were scanned using the NA_eff of 0.8 to achieve the higher resolution.

Fig. 1 k-space sampling trajectories. (a) Spiral scanning trajectories in the k-space. Red and green dots show the low and high NA_eff cases for 300 illumination angles, respectively. (b-c) Reconstructed 3-D k-space scattering potentials for (b) low (0.5) and (c) high (0.8) NA_eff, respectively.

Download Full Size | PDF

3.2. Reconstruction setup

We performed reconstruction using custom-made scripts in MatLab R2014a (MathWorks Inc., Natick, MA, USA) on a desktop computer (Intel Core i7-4770 CPU, 3.4 GHz, 12 GB RAM). To accelerate the ODT reconstruction speed, we utilized a graphics processing unit (GPU, Geforce GTX 660Ti, nVidia Corp., Santa Clara, CA, USA); custom-made functions based on the Compute Unified Device Architecture (CUDA) were used. Table 1 provides the information about the data size, parameters, and the computation time. Parameters were manually selected by trials and errors. The computation time depends on the number of interactions. In this study, the number of total iterations were 20 for the GP algorithm, and 100 for EP regularization (inner iteration: 10, outer iteration: 10) and the TV regularization (inner iteration: 20, outer iteration: 5). For a quantitative evaluation, we calculated the mean square error (MSE) of reconstruction images, defined as

M S E (x_{r e c o n}) = \frac{{‖ x_{t r u e} - x_{r e c o n} ‖}^{2}}{{‖ x_{t r u e} ‖}^{2}} .

Table 1. Description of data, used parameters, and the computation time

View Table

4. Results

4.1. Numerical simulation

We first performed numerical simulations with a phantom, a microsphere with a homogeneous RI of 1.44 and the diameter of 5 μm. We compared various reconstruction algorithms under two cases, NA_eff = 0.5 and 0.8, respectively. Figures 2(a)–2(d) show cross-sectional slices of the RI tomograms reconstructed with each algorithm (FT stands for the inverse 3-D Fourier transform of 3-D k-space filled according to diffraction algorithm with zero-padding). In the case of GP, the shape of a reconstructed RI tomogram along the z-axis is not confined even for the high NA_eff. For a more precise analysis, we considered the histogram of RI values and the RI profiles along the optical axis as shown in Figs. 2(e)–2(h). Although all three reconstruction algorithms show the increase in the RI value, the most frequent RI values are different and were obtained as 1.416 (GP), 1.43 (EP) and 1.44 (TV) as in Fig. 2(e) for the case of NA_eff = 0.5. In the case of NA_eff = 0.8, all three reconstruction algorithms generate RI tomograms with the RI peaks near 1.44 as in Fig. 2(f). While GP shows the increased RI values compared to FT, the distributions of reconstructed RI values are in broad ranges even for the case with high NA_eff. In contrast, TV shows the narrowest RI distributions. EP also shows more concentrated distributions compared to GP. Figures 2(g)–2(h) show the RI profiles along the optical axis. The width of RI profiles are 8.52 μm (GP), 6.80 μm (EP) and 6.49 μm (TV), respectively, as in Fig. 2(g) for the case of NA_eff = 0.5. In case of NA_eff = 0.8, the improved reduction of elongations are observed along the optical axis converging to the original value of RI (gold line) as in Fig. 2(h). GP shows the underestimation at the rim of the phantom. However, EP and TV show piecewise constant regions since they suppress noise preserving edges at the same time. Although EP results in the lowest MSE for the case of the high NA_eff, the results are dependent upon the NA_eff, which is not the case of TV.

Fig. 2 Numerical simulation results of a 5-μm-diameter fused silica bead phantom: 3-D RI tomograms numerically reconstructed by (a) FT, (b) GP, (c) EP, and (d) TV regularizations. (a-d); left column: z = 0 slice, right column: x = 0 slice, upper row: NA_eff = 0.5, lower row: NA_eff = 0.8, and the boundary of ROI is shown in white line. (e-f) Histograms of RI values in ROI divided by the total number of pixels in the ROI. (g-h) RI profiles along the z-axis. For figures (e-h), NA_eff = 0.5 [(e) and (g)] and NA_eff = 0.8 [(f) and (h)] were used. All scale bars are 2 μm.

Download Full Size | PDF

As a reference, we performed the simulation in case of noisy data in two cases (SNR = 5dB and 10dB, respectively). They are attached as the supplementary figure in Appendix B.

4.2. Fused silica microsphere

The experimental data of a fused silica microsphere were measured under two different NA_eff s, NA_eff = 0.5 and 0.8, respectively. Different slice images (xy and yz slices) of reconstructed RI tomograms for each algorithm are shown in Figs. 3(a–d). Reconstruction results using GP show the elongation of RI tomograms along the optical axis. Therefore, even though the reconstructed RI tomogram is reliable at the z = 0 plane, it shows defocused structures at high z planes. In particular, in Fig. 3(b), coherent optical noises are noticeable, which may cause difficulties in observing the intra-structures of a sample. The peak points of the histogram of RI distributions are 1.41 (GP), 1.412 (EP), and 1.44 (TV) as in Fig. 3(e) for the case of NA_eff = 0.5. With NA_eff = 0.8, RI values increase forming peak points at 1.438 (GP), 1.436 (EP) and 1.442 (TV) as shown in Fig. 3(f). Overall, reconstructed RI values are lower compared to the previous numerical simulation. Also, GP and EP do not converge even for the high NA_eff case; in the case of GP, the distributions of RI tomograms become wider with the increase of the NA_eff. We also compared the width of reconstructed RI profiles along the optical axis in Fig. 3(g) and 3(h). As the NA_eff increases, the width changes from 8.62 μm to 6.20 μm (GP), from 8.42 μm to 5.48 μm (EP), and from 6.29 μm to 4.97 μm (TV), respectively. As shown in the numerical simulation, EP and TV regularization show less fluctuations in the RI tomograms compared to GP. TV shows the most piecewise constant regions and the lowest dependency on the NA_eff.

Fig. 3 Experimental results of a 5-μm-diameter fused silica bead: 3-D RI tomograms numerically reconstructed by (a) FT, (b) GP, (c) EP, and (d) TV regularizations. (a-d); left column: z = 0 slice, right column: x = 0 slice, upper row: NA_eff = 0.5, lower row: NA_eff = 0.8, and the boundary of ROI is shown in white line. (e-f) Histograms of RI values in ROI divided by the total number of pixels in the ROI. (g-h) RI profiles along the z-axis. For figures (e-h), NA_eff = 0.5 [(e) and (g)] and NA_eff = 0.8 [(f) and (h)] were used. All scale bars are 2 μm.

Download Full Size | PDF

4.3. Biological samples:RBCs and hepatocytes

Thanks to its capability of 3-D stainless visualization, the ODT has been utilized for imaging 3-D structures of live biological cells, and quantification of the physical characteristics of the individual cells such as shape, volume, and mass [1, 2, 61]. Therefore, in this subsection, we compared the algorithms using real biological cells: RBCs and hepatocytes. Figures 4 and 5 illustrate the reconstruction results. RBCs should have narrow RI distribution, since RBC does not have subcellular structures. However, the central dimple region of a RBC reconstruction is thinner due to the loss of information at the missing cone, so the reconstructed distribution is not uniform [40]. This problem can be remedied with a higher NA_eff such as a water-immersed objective lens [62].

Fig. 4 Experimental results of RBC RI tomograms obtained with (a) FT, (b) GP, (c) EP, and (d) TV regularization. The white dotted lines represent the slices of the complementary figures. All scale bars are 3 μm.

Download Full Size | PDF

Fig. 5 Experimental results of hepatocyte RI tomograms obtained with (a) FT, (b) GP, (c) EP, and (d) TV regularization. For figures (a-d), the white dotted lines represent the slices of the complementary figures. All scale bars are 5 μm.

Download Full Size | PDF

For real biological samples, the noise-amplifying drawbacks of the GP algorithm can be more significant due to movements, fluctuations, and floating debris of living cells as in Figs. 4(b) and 5(b). The EP shows the capacity of denoising and preserving the texture at the same time, as observed in the texture of RBC in Fig. 4(c) and the cytoplasm region of hepatocyte in Fig. 5(c). TV regularized reconstruction shows the highest RI value similar to the previous microsphere case. However, in Fig. 4(d) and Fig. 5(d), TV results show the cartoon-like artifacts. Cartoon-like artifacts appear when we impose the total variation regularization on the cost function too much. It means that an image becomes cartoon-like with sharp jumps between nearly piecewise constant regions which appears in Fig. 4(d) [63]. Especially in case of Fig. 5(d), we can see the stair-shaped artifact which appears in the cytoplasm region of the hepatocyte. Therefore, we found that TV regularization is not appropriate for samples with inhomogeneous RI distributions, such as cytoplasms of hepatocyte cells. EP results, however, show reasonable RI values and textures compared to the raw tomogram.

5. Conclusion

In this paper, we conducted the systematic comparative study on the effects of existing regularization algorithms to 3-D RI tomograms. GP, EP, and TV algorithms were compared both numerically and experimentally using 5-μm-diameter fused silica microspheres, RBCs, and hepatocytes. Using GP, EP, and TV, RI values increase. However, GP shows the limitation in the mitigation of the elongation of RI tomograms and noise handling capacity. The TV regularization shows excellent reconstruction performance for samples with relatively homogenous RI values, even for low NA_eff condition. However, RI tomograms with many textures could not be reconstructed correctly using TV regularization. The EP regularization is more efficient than TV even though EP regularization shows higher dependency of the NA_eff.

Appendix A

We can rewrite the equation Eq. (14) by splitting variables

(f^{m + 1}, d^{m + 1}, v^{m + 1}) = \min_{f, d, v} | d | + I (v \geq 0) + \frac{μ}{2} {‖ A f - g^{k} ‖}_{2}^{2} such that d = d (f), v = f .

The associated unconstrained cost function is

(f^{m + 1}, d^{m + 1}, v^{m + 1}) = \min_{f, d, v} | d | + I (v \geq 0) + \frac{μ}{2} {‖ A f - g^{k} ‖}_{2}^{2} + \frac{α}{2} ‖ d - D (f) ‖_{2}^{2} + \frac{β}{2} ‖ v - f ‖_{2}^{2} .

By specifying d_x = ∇_x f, d_y = ∇_y f and d_z = ∇_z f, the split Bregman iterations are given by,

\begin{array}{l} (f^{m + 1}, d_{x}^{m + 1}, d_{y}^{m + 1}, d_{z}^{m + 1}, v^{m + 1}) \\ = \min_{f, d, v} {‖ d_{x}, d_{y}, d_{z} ‖}_{2} + I (v \geq 0) + \frac{μ}{2} {‖ A f - g^{k} ‖}_{2}^{2} + \frac{α}{2} ‖ d_{x} - \nabla_{x} f - b_{x}^{m} ‖_{2}^{2} \\ + \frac{α}{2} ‖ d_{y} - \nabla_{y} f - b_{y}^{m} ‖_{2}^{2} + \frac{α}{2} ‖ d_{z} - \nabla_{z} f - b_{z}^{m} ‖_{2}^{2} + \frac{β}{2} ‖ v - f - b_{v}^{m} ‖_{2}^{2}, \\ b_{x}^{m + 1} = b_{x}^{m} + (\nabla_{x} (f^{m + 1}) - d_{x}^{m + 1}), \\ b_{y}^{m + 1} = b_{y}^{m} + (\nabla_{y} (f^{m + 1}) - d_{y}^{m + 1}), \\ b_{z}^{m + 1} = b_{z}^{m} + (\nabla_{z} (f^{m + 1}) - d_{z}^{m + 1}), \\ b_{v}^{m + 1} = b_{v}^{m} + (f^{m + 1}) - v^{m + 1}), \end{array}

where

{‖ d_{x}, d_{y}, d_{z} ‖}_{2} = \sum_{n} \sqrt{d_{x, n}^{2} + d_{y, n}^{2} + d_{z, n}^{2}}

. To calculate f^m+¹, we have to solve the subproblem

\begin{matrix} \min_{f} \frac{μ}{2} {‖ A f - g^{k} ‖}_{2}^{2} + \frac{α}{2} ‖ d_{x}^{m} - \nabla_{x} f - b_{x}^{m} ‖_{2}^{2} + \frac{α}{2} ‖ d_{y}^{m} - \nabla_{y} f - b_{y}^{m} ‖_{2}^{2} \\ + \frac{α}{2} ‖ d_{z}^{m} - \nabla_{z} f - b_{z}^{m} ‖_{2}^{2} + \frac{β}{2} ‖ v_{m} - f - b_{v}^{m} ‖_{2}^{2} . \end{matrix}

By differentiating it, we can get the following equation

(μ A^{H} A - α Δ + β I) f^{m}^{+ 1} = r h s^{m},

where

Δ = - (\nabla_{x}^{T} \nabla_{x} + \nabla_{y}^{T} \nabla_{y} + \nabla_{z}^{T} \nabla_{z})

and

r h s^{m} = μ A^{H} g^{k} + α \nabla_{x}^{T} (d_{x}^{m} - b_{x}^{m}) + α \nabla_{y}^{T} (d_{y}^{m} - b_{y}^{m}) + α \nabla_{z}^{T} (d_{z}^{m} - b_{z}^{m}) + β (v^{m} - b_{v}^{m})

. To calculate d_x,d_y and d_z, we have to solve subproblems

\min_{d} ‖ d_{x}, d_{y}, d_{z} |_{2} + \frac{α}{2} ‖ d_{x} - \nabla_{x} f^{m + 1} - b_{x}^{m} ‖_{2}^{2} + \frac{α}{2} ‖ d_{y} - \nabla_{y} f^{m + 1} - b_{y}^{m} ‖_{2}^{2} + \frac{α}{2} ‖ d_{z} - \nabla_{z} f^{m + 1} - b_{z}^{m} ‖_{2}^{2} .

The equation Eq. (22) can be minimized with respect to d_x,d_y and d_z respectively using a generalized shrinkage formula.

\begin{matrix} d_{x}^{m + 1} = \max (s^{m} - \frac{1}{α}, 0) \frac{\nabla_{x} f^{m + 1} + b_{x}^{m}}{s^{m}}, \\ d_{y}^{m + 1} = \max (s^{m} - \frac{1}{α}, 0) \frac{\nabla_{y} f^{m + 1} + b_{y}^{m}}{s^{m}}, \\ d_{z}^{m + 1} = \max (s^{m} - \frac{1}{α}, 0) \frac{\nabla_{z} f^{m + 1} + b_{z}^{m}}{s^{m}}, \end{matrix}

where

s^{m} = \sqrt{| \nabla_{x} f^{m}^{+ 1} + b_{x}^{m} |^{2} + | \nabla_{y} f^{m}^{+ 1} + b_{y}^{m} |^{2} + | \nabla_{z} f^{m}^{+ 1} + b_{z}^{m} |^{2}}

. Finally, we can get v^m⁺¹ by solving the subproblem

\min_{v} I (v \geq 0) + \frac{β}{2} ‖ v - f^{m + 1} - b_{v}^{m} ‖_{2}^{2},

which corresponds to the proximal mapping of indicator function [54]

v^{m + 1} = \max (f^{m + 1} + b_{v}^{m}, 0) .

Appendix B

We added noise in fused silica microsphere simulation data with the different levels of SNR: 5dB and 10dB, as shown in Fig. 6 and Fig. 7, respectively. The overall tendency is similar with the noiseless case. Although GP algorithm shows the limitation in noise handling capacity considering MSE, EP and TV show better results in denoising.

Fig. 6 Numerical simulation results of a 5-μm-diameter fused silica bead phantom with noise (SNR = 5dB): 3-D RI tomograms numerically reconstructed by (a) FT, (b) GP, (c) EP, and (d) TV regularizations. (a-d); left column: z = 0 slice, right column: x = 0 slice, upper row: NA_eff = 0.5, lower row: NA_eff = 0.8, and the boundary of ROI is shown in white line.

Download Full Size | PDF

Fig. 7 Numerical simulation results of a 5-μm-diameter fused silica bead phantom with noise (SNR = 10dB): 3-D RI tomograms numerically reconstructed by (a) FT, (b) GP, (c) EP, and (d) TV regularizations. (a-d); left column: z = 0 slice, right column: x = 0 slice, upper row: NA_eff = 0.5, lower row: NA_eff = 0.8, and the boundary of ROI is shown in white line.

Download Full Size | PDF

Acknowledgments

This work was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIP) (No. NRF-2014R1A2A1A11052491).

References and links

1. G. Popescu, Quantitative Phase Imaging of Cells and Tissues (McGraw Hill Professional, 2011).

2. K. Lee, K. Kim, J. Jung, J. Heo, S. Cho, S. Lee, G. Chang, Y. Jo, H. Park, and Y. Park, “Quantitative phase imaging techniques for the study of cell pathophysiology: from principles to applications,” Sensors 13, 4170–4191 (2013). [CrossRef] [PubMed]

3. G. Popescu, T. Ikeda, R. R. Dasari, and M. S. Feld, “Diffraction phase microscopy for quantifying cell structure and dynamics,” Opt. Lett. 31, 775–777 (2006). [CrossRef] [PubMed]

4. T. Ikeda, G. Popescu, R. R. Dasari, and M. S. Feld, “Hilbert phase microscopy for investigating fast dynamics in transparent systems,” Opt. Lett. 30, 1165–1167 (2005). [CrossRef] [PubMed]

5. Y. Park, G. Popescu, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Diffraction phase and fluorescence microscopy,” Opt. Express 14, 8263–8268 (2006). [CrossRef] [PubMed]

6. Y. Park, W. Choi, Z. Yaqoob, R. Dasari, K. Badizadegan, and M. S. Feld, “Speckle-field digital holographic microscopy,” Opt. Express 17, 12285–12292 (2009). [CrossRef] [PubMed]

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

8. B. Simon, M. Debailleul, A. Beghin, Y. Tourneur, and O. Haeberlé, “High-resolution tomographic diffractive microscopy of biological samples,” J. Biophotonics 3, 462–467 (2010). [CrossRef] [PubMed]

9. K. Kim, H. Yoon, M. Diez-Silva, M. Dao, R. R. Dasari, and Y. Park, “High-resolution three-dimensional imaging of red blood cells parasitized by plasmodium falciparum and in situ hemozoin crystals using optical diffraction tomography,” J. Biomed. Opt. 19, 011005 (2014). [CrossRef]

10. Y. Park, C. A. Best, T. Auth, N. S. Gov, S. A. Safran, G. Popescu, S. Suresh, and M. S. Feld, “Metabolic remodeling of the human red blood cell membrane,” Proc. Natl. Acad. Sci. U.S.A. 107, 1289–1294 (2010). [CrossRef] [PubMed]

11. Y. Park, C. A. Best, K. Badizadegan, R. R. Dasari, M. S. Feld, T. Kuriabova, M. L. Henle, A. J. Levine, and G. Popescu, “Measurement of red blood cell mechanics during morphological changes,” Proc. Natl. Acad. Sci. USA 107, 6731–6736 (2010). [CrossRef] [PubMed]

12. Y. Park, C. A. Best, T. Kuriabova, M. L. Henle, M. S. Feld, A. J. Levine, and G. Popescu, “Measurement of the nonlinear elasticity of red blood cell membranes,” Phys. Rev. E 83, 051925 (2011). [CrossRef]

13. G. Popescu, Y. Park, N. Lue, C. Best-Popescu, L. Deflores, R. R. Dasari, M. S. Feld, and K. Badizadegan, “Optical imaging of cell mass and growth dynamics,” Am. J. Physiol. Cell Physiol. 295, C538–C544 (2008). [CrossRef] [PubMed]

14. M. Mir, Z. Wang, Z. Shen, M. Bednarz, R. Bashir, I. Golding, S. G. Prasanth, and G. Popescu, “Optical measurement of cycle-dependent cell growth,” Proc. Natl. Acad. Sci. USA 108, 13124–13129 (2011). [CrossRef] [PubMed]

15. Y. Park, T. Yamauchi, W. Choi, R. Dasari, and M. S. Feld, “Spectroscopic phase microscopy for quantifying hemoglobin concentrations in intact red blood cells,” Opt. Lett. 34, 3668–3670 (2009). [CrossRef] [PubMed]

16. Y. Jang, J. Jang, and Y. Park, “Dynamic spectroscopic phase microscopy for quantifying hemoglobin concentration and dynamic membrane fluctuation in red blood cells,” Opt. Express 20, 9673–9681 (2012). [CrossRef] [PubMed]

17. J. Jung, K. Kim, H. Yu, K. Lee, S. Lee, S. Nahm, H. Park, and Y. Park, “Biomedical applications of holographic microspectroscopy [invited],” Appl. Opt. 53, G111–G122 (2014). [CrossRef] [PubMed]

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

19. M. Fauver, E. Seibel, J. R. Rahn, M. Meyer, F. Patten, T. Neumann, and A. Nelson, “Three-dimensional imaging of single isolated cell nuclei using optical projection tomography,” Opt. Express 13, 4210–4223 (2005). [CrossRef] [PubMed]

20. G. G. Levin, G. N. Vishnyakov, C. S. Zakarian, A. V. Likhachov, V. V. Pickalov, G. I. Kozinets, J. K. Novoderzhkina, and E. A. Streletskaya, “Three-dimensional limited-angle microtomography of blood cells: experimental results,” in “BiOS’98 International Biomedical Optics Symposium,” (International Society for Optics and Photonics, 1998), pp. 159–164.

21. V. Lauer, “New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope,” J. Microsc. 205, 165–176 (2002). [CrossRef] [PubMed]

22. B. Simon, M. Debailleul, V. Georges, V. Lauer, and O. Haeberlé, “Tomographic diffractive microscopy of transparent samples,” Eur. Phys. J. Appl. Phys. 44, 29–35 (2008). [CrossRef]

23. S. O. Isikman, W. Bishara, S. Mavandadi, W. Y. Frank, S. Feng, R. Lau, and A. Ozcan, “Lens-free optical tomographic microscope with a large imaging volume on a chip,” Proc. Natl. Acad. Sci. USA 108, 7296–7301 (2011). [CrossRef] [PubMed]

24. L. Yu and M. K. Kim, “Wavelength-scanning digital interference holography for tomographic three-dimensional imaging by use of the angular spectrum method,” Opt. Lett. 30, 2092–2094 (2005). [CrossRef] [PubMed]

25. J. Kühn, F. Montfort, T. Colomb, B. Rappaz, C. Moratal, N. Pavillon, P. Marquet, and C. Depeursinge, “Submicrometer tomography of cells by multiple-wavelength digital holographic microscopy in reflection,” Opt. Lett. 34, 653–655 (2009). [CrossRef] [PubMed]

26. M. Potcoava and M. Kim, “Optical tomography for biomedical applications by digital interference holography,” Meas. Sci. Technol. 19, 074010 (2008). [CrossRef]

27. J. A. Izatt, E. A. Swanson, J. G. Fujimoto, M. R. Hee, and G. M. Owen, “Optical coherence microscopy in scattering media,” Opt. Lett. 19, 590–592 (1994). [CrossRef] [PubMed]

28. M. A. Choma, A. K. Ellerbee, C. Yang, T. L. Creazzo, and J. A. Izatt, “Spectral-domain phase microscopy,” Opt. Lett. 30, 1162–1164 (2005). [CrossRef] [PubMed]

29. Z. Wang, D. L. Marks, P. S. Carney, L. J. Millet, M. U. Gillette, A. Mihi, P. V. Braun, Z. Shen, S. G. Prasanth, and G. Popescu, “Spatial light interference tomography (slit),” Opt. Express 19, 19907–19918 (2011). [CrossRef] [PubMed]

30. O. Haeberlé, K. Belkebir, H. Giovaninni, and A. Sentenac, “Tomographic diffractive microscopy: basics, techniques and perspectives,” J. Mod. Opt. 57, 686–699 (2010). [CrossRef]

31. R. Fiolka, K. Wicker, R. Heintzmann, and A. Stemmer, “Simplified approach to diffraction tomography in optical microscopy,” Opt. Express 17, 12407–12417 (2009). [CrossRef] [PubMed]

32. S. S. Kou and C. J. Sheppard, “Image formation in holographic tomography: high-aperture imaging conditions,” Appl. Opt. 48, H168–H175 (2009). [CrossRef] [PubMed]

33. Y. Ruan, P. Bon, E. Mudry, G. Maire, P. Chaumet, H. Giovannini, K. Belkebir, A. Talneau, B. Wattellier, S. Mon-neret, and A. Sentenac, “Tomographic diffractive microscopy with a wavefront sensor,” Opt. Lett. 37, 1631–1633 (2012). [CrossRef] [PubMed]

34. S. Uttam, S. A. Alexandrov, R. K. Bista, and Y. Liu, “Tomographic imaging via spectral encoding of spatial frequency,” Opt. Express 21, 7488–7504 (2013). [CrossRef] [PubMed]

35. A. Kuś, M. Dudek, B. Kemper, M. Kujawińska, and A. Vollmer, “Tomographic phase microscopy of living three-dimensional cell cultures,” J. Biomed. Opt. 19, 046009 (2014). [CrossRef]

36. L. Su, L. Ma, and H. Wang, “Improved regularization reconstruction from sparse angle data in optical diffraction tomography,” Appl. Opt. 54, 859–868 (2015). [CrossRef] [PubMed]

37. L. Ma, H. Wang, L. Su, Y. Li, and H. Jin, “Digital holographic microtomography with few angle data-sets,” J. Mod. Opt. 61, 1140–1146 (2014). [CrossRef]

38. F. Merola, L. Miccio, P. Memmolo, G. Di Caprio, A. Galli, R. Puglisi, D. Balduzzi, G. Coppola, P. Netti, and P. Ferraro, “Digital holography as a method for 3d imaging and estimating the biovolume of motile cells,” Lab Chip 13, 4512–4516 (2013). [CrossRef] [PubMed]

39. Y. Cotte, F. Toy, P. Jourdain, N. Pavillon, D. Boss, P. Magistretti, P. Marquet, and C. Depeursinge, “Marker-free phase nanoscopy,” Nature Photon. 7, 113–117 (2013). [CrossRef]

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

41. Y.-C. Lin and C.-J. Cheng, “Sectional imaging of spatially refractive index distribution using coaxial rotation digital holographic microtomography,” J. Opt. 16, 065401 (2014). [CrossRef]

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

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

44. S. Bartolac, R. Clackdoyle, F. Noo, J. Siewerdsen, D. Moseley, and D. Jaffray, “A local shift-variant fourier model and experimental validation of circular cone-beam computed tomography artifacts,” Med. Phys. 36, 500–512 (2009). [CrossRef] [PubMed]

45. M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magn. Reson. Med. 58, 1182–1195 (2007). [CrossRef] [PubMed]

46. R. Gerchberg, “Super-resolution through error energy reduction,” J. Mod. Opt. 21, 709–720 (1974).

47. A. Papoulis, “A new algorithm in spectral analysis and band-limited extrapolation,” IEEE Trans. Circuits and Systems 22, 735–742 (1975). [CrossRef]

48. P. Charbonnier, L. Blanc-Féraud, G. Aubert, and M. Barlaud, “Deterministic edge-preserving regularization in computed imaging,” IEEE Trans. on Image Processing 6, 298–311 (1997). [CrossRef]

49. A. H. Delaney and Y. Bresler, “Globally convergent edge-preserving regularized reconstruction: an application to limited-angle tomography,” IEEE Trans. Image Processing 7, 204–221 (1998). [CrossRef]

50. E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Information Theory 52, 489–509 (2006). [CrossRef]

51. Y. Sung, W. Choi, N. Lue, R. R. Dasari, and Z. Yaqoob, “Stain-free quantification of chromosomes in live cells using regularized tomographic phase microscopy,” PLoS ONE 7, e49502 (2012). [CrossRef] [PubMed]

52. S. J. LaRoque, E. Y. Sidky, and X. Pan, “Accurate image reconstruction from few-view and limited-angle data in diffraction tomography,” J. Opt. Soc. Am. A 25, 1772–1782 (2008). [CrossRef]

53. T. Goldstein and S. Osher, “The split Bregman method for l1-regularized problems,” SIAM. J. Imaging. Sci. 2, 323–343 (2009). [CrossRef]

54. J.-J. Abascal, J. Chamorro-Servent, J. Aguirre, S. Arridge, T. Correia, J. Ripoll, J. J. Vaquero, and M. Desco, “Fluorescence diffuse optical tomography using the split Bregman method,” Med. Phys. 38, 6275–6284 (2011). [CrossRef] [PubMed]

55. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (Society for Industrial and Applied Mathematics, 2001). [CrossRef]

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

57. L. M. Bregman, “The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming,” USSR. Comp. Math. Math+. 7, 200–217 (1967). [CrossRef]

58. Y. Kim, H. Shim, K. Kim, H. Park, J. H. Heo, J. Yoon, C. Choi, S. Jang, and Y. Park, “Common-path diffraction optical tomography for investigation of three-dimensional structures and dynamics of biological cells,” Opt. Express 22, 10398–10407 (2014). [CrossRef] [PubMed]

59. O. Zhernovaya, O. Sydoruk, V. Tuchin, and A. Douplik, “The refractive index of human hemoglobin in the visible range,” Phys. Med. Biol. 56, 4013 (2011). [CrossRef] [PubMed]

60. K. Kim, K. S. Kim, H. Park, J. C. Ye, and Y. Park, “Real-time visualization of 3-D dynamic microscopic objects using optical diffraction tomography,” Opt. Express 21, 32269–32278 (2013). [CrossRef]

61. Y. Kim, H. Shim, K. Kim, H. Park, S. Jang, and Y. Park, “Profiling individual human red blood cells using common-path diffraction optical tomography,” Sci. Rep. 4, 6659 (2014). [CrossRef] [PubMed]

62. K. Kim, Z. Yaqoob, K. Lee, J. W. Kang, Y. Choi, P. Hosseini, P. T. So, and Y. Park, “Diffraction optical tomography using a quantitative phase imaging unit,” Opt. Lett. 39, 6935–6938 (2014). [CrossRef] [PubMed]

63. P. Getreuer, “Rudin-Osher-Fatemi total variation denoising using split Bregman,” Image Processing on Line 10, 74 (2012).

Data (single type) Size(x×y×z)	NA_eff	GP	EP				TV
Data (single type) Size(x×y×z)	NA_eff	Time (s)	α	β	T	Time (s)	μ	α	β	Time (s)
phantom 343 × 320 × 170	NA_eff = 0.5	15.80	1	1	0.2	170.79	2	2	0.1	408.60
phantom 343 × 320 × 170	NA_eff = 0.8	16.88	1	1	0.12	170.43	0.6	0.6	0.1	407.38
microsphere 343 × 320 × 170	NA_eff = 0.5	16.47	4	5	0.1	170.37	3	3	0.5	408.29
microsphere 343 × 320 × 170	NA_eff = 0.8	17.83	2	3	0.3	169.48	0.5	0.5	0.5	409.15
RBC 256 × 256 × 256 hepatocyte	NA_eff = 0.8	12.93	1	1	0.2	145.59	1	1	10	364.77
374 × 350 × 170	NA_eff = 0.8	22.27	1.5	1	0.15	221.60	1	1	1	506.65

Comparative study of iterative reconstruction algorithms for missing cone problems in optical diffraction tomography

Abstract

1. Introduction

2. Theory

2.1. Optical diffraction tomography

2.2. Problem formulation

2.3. Iterative reconstruction methods

2.3.1. Convergent EP regularization with the non-negativity constraint

2.3.2. TV regularization with non-negativity constraint

3. Method

3.1. Experimental setup

3.2. Reconstruction setup

4. Results

4.1. Numerical simulation

4.2. Fused silica microsphere

4.3. Biological samples:RBCs and hepatocytes

5. Conclusion

Appendix A

Appendix B

Acknowledgments

References and links

Cited By

Figures (7)

Tables (1)

Equations (25)

Optics Express