Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Iterative optical diffraction tomography for illumination scanning configuration

Open Access Open Access

Abstract

Optical diffraction tomography (ODT) is used to reconstruct refractive-index distributions from multiple measurements in the object rotating configuration (ORC) or the illumination scanning configuration (ISC). Because of its fast data acquisition and stability, ISC-based ODT has been widely used for biological imaging. ODT typically fails to reconstruct multiply-scattering samples. The previously developed iterative ODT (iODT) was for the multiply-scattering objects in ORC, and could not be directly applied to ISC. To resolve this mismatch, we developed an ISC update and numerically demonstrated its accuracy. With the same prior knowledge, iODT-ISC outperforms conventional ODT in resolving the missing-angle problem.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

As a label-free and non-invasive imaging technique, optical diffraction tomography (ODT), has been widely used in imaging biological samples with sub-wavelength resolution [19]. Conventional ODT reconstructs the refractive index (RI) distribution of the sample based on the diffracted fields over multiple projection views. The diffracted fields are holographically measured by use of the phase-shifting or off-axis interferometry [10,11]. To obtain the fields diffracted by the sample over many angles, two configurations are used: the object rotating configuration (ORC) [1216] and the illumination scanning configuration (ISC) [1,8,9].

Originally, ODT was developed for ORC where the sample is rotated, while the probing wave and camera are fixed. Since this configuration is relatively simple and inexpensive, it has been widely employed in imaging stationary samples, such as optical fibers [12]. The sample is typically rotated by a motorized rotation stage, which may introduce mechanical instability. To apply ORC for bio-imaging, several alternative rotation configurations were recently invented based on optical tweezers [17], dielectrophoretic cell rotation [18], natural flow of cells [19], and acoustic pressure [20]; however, the complexity or cost of the setup can be increased. By contrast, the ISC rotates the probing wave by use of a beam steering device (e.g., galvanometer, digital micromirror device, or spatial light modulator) [1,8,21,22], while keeping the sample and the camera fixed. ODT based on ISC has good transverse resolution, but poor axial resolution caused by inadequate filling of the spectrum along the axial direction [23]. Because of its stability and high-speed data acquisition, ISC has been frequently used in commercial bio-imaging instruments (e.g., by Nanolive, ltd., and Tomocube, Inc.).

Conventional ODT relies on the Fourier diffraction theorem, which is applicable under the condition of weak scattering [24]. However, many biological samples, such as clusters of cancer cells, have high contrasts and large optical path-length differences (OPDs) and therefore tend to produce multiple scattering. When the object produces multiple scattering or is surrounded by multiply-scattering media, ODT becomes ineffective and the inversion problem becomes more challenging. To extend ODT application to multiply-scattering samples, optimization-based ODT methods were developed in which multiple scattering can be modelled by beam propagation method [9,25,26], wave propagation method [27,28] or Lippmann-Schwinger equation [29,30]; however, the optimization methods can be time consuming and may be trapped in local minima, especially for objects introducing large OPDs. More recently, a more computationally efficient iterative optical diffraction tomography (iODT) was proposed for the reconstruction beyond the weakly-scattering regime [31,32].

The frequency domain mappings for ORC and ISC are based on the rotation and shifting of the Ewald spheres, respectively. This results in different shapes and densities of the frequency supports, implying that the over-representation of the Ewald spheres at the same frequency component needs to be handled differently. Therefore, iODT as originally developed for ORC, has artifacts if directly applied for the reconstruction from ISC measurements. Therefore, a method to extend iODT for ISC is desired. In this paper, we propose a new iODT algorithm for ISC-based measurements (hereafter referred to as iODT-ISC) that deals with the configuration mismatch.

In iODT-ISC reconstruction, the incident and measured fields are respectively forward-backward propagated through an estimate of the object for each of the titled illumination angles. The perturbative Rytov phases computed from the resultant propagated fields are summed in the spatial domain and normalized in the frequency domain by the appropriate coherent transfer function, to form a perturbative correction to the current estimate of the object. Numerical results demonstrate that the proposed iODT-ISC technique provides accurate and efficient reconstructions of multiply-scattering objects from ISC-based measurements.

2. Theory

For simplicity, two-dimensional (2D) samples, such as cross sections of optical fibers, are considered in this paper. The sample to be imaged has a RI distribution $n(\overrightarrow r ) = {n_\textrm{b}} + \Delta n(\overrightarrow r )$, where ${n_\textrm{b}}$ is the RI of the background medium and $\overrightarrow r = (x,y)$ is the spatial coordinate. The object function (viz. scattering potential) is proportional to the permittivity difference between the sample and background medium, i.e., $f(\overrightarrow r ) = {(2\pi {k_0})^2} \cdot [{n^2}(\overrightarrow r ) - n_\textrm{b}^2]$, where ${k_0} = 1/{\lambda _0}$ is the wavenumber, ${\lambda _0}$ is the wavelength of the probing wave in free-space. A probing wave passing through the sample produces a diffracted field which is collected by the imaging system and interferometrically detected by a camera using the phase-shifting or off-axis holography. The diffracted fields corresponding to different angles of incidence on the sample form a sinogram $u(y;{\theta _{\textrm{ill}}})$.

The ORC and ISC configurations for ODT are illustrated in Fig. 1. This paper develops iODT reconstruction algorithm for ISC, which is more suitable for bioimaging. We first discuss two weakly-scattering reconstruction methods, respectively implemented in the frequency and spatial domains, and then introduce iODT-ISC reconstruction, which is applicable to multiply-scattering objects.

 figure: Fig. 1.

Fig. 1. Configurations for ODT (a) object rotating configuration (ORC) and (b) illumination scanning configuration (ISC).

Download Full Size | PDF

2.1 Direct Fourier interpolation (DFI)

Figure 2 illustrates the DFI method, which is implemented in the frequency domain. Since the relay imaging system conjugates the plane of the camera to the plane of best focus in the object, the sinogram interferometrically measured by the camera is the focused sinogram. We employ the Rytov approximation which has been shown to be superior to the Born approximation [33]. For a given illumination angle, a sliced spectrum of the object function is proportional to the Fourier transform of the quantity related to the refocused diffracted field by the Fourier diffraction theorem,

$$\widehat F(\overrightarrow k - {\overrightarrow k _{\textrm{inc}}};{\theta _{\textrm{ill}}}) ={-} j4\pi {k_x} \cdot {\textbf{\textit{F}}}\left( {{u_{\textrm{bgd}}}(y;{\theta_{\textrm{ill}}})\ln \frac{{u(y;{\theta_{\textrm{ill}}})}}{{{u_{\textrm{bgd}}}(y;{\theta_{\textrm{ill}}})}}} \right),$$
where $\widehat F$ is the spectrum of the reconstructed object function, ${\overrightarrow k _{\textrm{inc}}}$ and $\overrightarrow k $ are the wavevectors of the incident and scattered fields, respectively, ${k_x} = \sqrt {k_\textrm{b}^2 - k_y^2} $ is the axial component of the scattered wavevector, ${\textbf{\textit{F}}}$ is the Fourier transform operator, $u(y;{\theta _{\textrm{ill}}})$ and ${u_{\textrm{bgd}}}(y;{\theta _{\textrm{ill}}})$ are the refocused diffracted fields and background fields, respectively. Once all the sliced spectra are mapped to the Ewald spheres specified by the illumination angles and averaged in the overlapped regions, the reconstructed object function can be obtained by use of the inverse Fourier transform.

 figure: Fig. 2.

Fig. 2. Illustration of direct Fourier interpolation (DFI) for ISC.

Download Full Size | PDF

2.2 Extended-depth-of-focus filter backpropagation for ISC (ISC-EDOF-FBPP)

ISC-FBPP is a spatial domain implementation of DFI, while ISC-EDOF-FBPP [34] introduces EDOF feature on top of ISC-FBPP. The over-representation of the same frequency components in the overlapping region of Ewald spheres is normalized by the coherent transfer function (CTF) [35]. The relay imaging system conjugates the camera plane with the plane of the best focus, and the depth of focus is very limited. To extend the depth of focus, the measured fields $u(x = 0,y;{\theta _{\textrm{ill}}})$ are forward-backward propagated, in the background medium along the illumination direction, to form the sinogram with EDOF $u(x,y;{\theta _{\textrm{ill}}})$. Then the complex Rytov phase

$${\phi _\textrm{R}}(\overrightarrow r ;{\theta _{\textrm{ill}}}) = \ln \frac{{u(\overrightarrow r ;{\theta _{\textrm{ill}}})}}{{{u_{\textrm{bgd}}}(\overrightarrow r ;{\theta _{\textrm{ill}}})}},$$
is calculated and summed in the special domain over all the illumination angles to form an “object-like” function
$$\widehat {f^{\prime}}(\overrightarrow r ) = \frac{1}{{ - j4\pi {k_\textrm{b}}}}\sum\nolimits_{{\theta _{\textrm{ill}}}} {{\phi _\textrm{R}}(\overrightarrow r ;{\theta _{\textrm{ill}}})\,\Delta {\theta _{\textrm{ill}}}} .$$
The function $\widehat {f^{\prime}}(\overrightarrow r )$ shows geometric features of the object approximately, such as the shape and location. The sum of Rytov phases in the spatial domain is associated with the sum of the Ewald spheres in the frequency domain, which results in an over-representation of the spectrum, especially in the vicinity of the low frequency component. Therefore, the “object-like” function $\widehat {f^{\prime}}(\overrightarrow r )$ needs to be normalized in the frequency domain by CTF, so as to obtain a more accurate reconstruction,
$$\widehat f(\overrightarrow r ) = {{\textbf{\textit{F}}}^{ - 1}}[{{\textbf{\textit{F}}}(\widehat {f^{\prime}}(\overrightarrow r ))/\textrm{CTF(}\overrightarrow k \textrm{)}} ].$$

2.3 Iterative optical diffraction tomography for ISC (iODT-ISC)

Since both DFI and ISC-EDOF-FBPP assume weak scattering, they are prone to failure when reconstructing an object with high contrast, large OPD, and complicated structure. Although our previously proposed iODT method can extend ODT to the multiply-scattering regime, it is not directly applicable to the ISC configuration. We have therefore developed a new algorithm iODT-ISC, summarized in Algorithm 1. The iterative process of iODT-ISC is similar to the ORC-based iODT algorithm, but it differs in two aspects: tilted beam propagation and computation of the perturbative correction of the object function from the field discrepancies. The tilted beam propagation method (BPM) simulates the optical field propagation slice by slice, in which the refraction and diffraction steps are modelled separately. The refraction step is formulated as a phase modulation ${\textrm{exp}} [j2\pi {k_0}\,\Delta n(y) \cdot \Delta x/\cos {\theta _{\textrm{ill}}}]$, where $\Delta n(y)$ and ${\theta _{\textrm{ill}}}$ are the RI variation at the current slice and the illumination angle, respectively. The tilted diffraction step is formulated by the modified shifted angular spectrum method [36] or the translated angular spectrum method [37]. The computation of the perturbative correction in iODT-ISC will be explained as follows, with the emphasis on its difference from that in the original iODT.

The inputs of the iODT-ISC algorithm are the known incident fields ${u_0}(x ={-} d,y;{\theta _{\textrm{ill}}})$ and the true diffracted fields ${u_\textrm{t}}(x = d,y;{\theta _{\textrm{ill}}})$ at the boundaries of simulation dimensions for each illumination angle ${\theta _{\textrm{ill}}}$, the initial guess of the object function ${f^{(0)}}$(typically zero), and the pre-defined maximum number of iterations ${m_{\max }}$. If prior information exists, then the constraint ${\textbf{\textit{C}}}$ can be applied to the reconstructed object function. The known incident fields and the true diffracted fields are respectively forward and backward propagated through the current estimate of the object function, to form the forward and backward fields for all the illumination angles. These forward and backward propagated fields produce perturbative Rytov phases. In the original iODT (for ORC), the perturbative Rytov phases are high-pass filtered in the frequency domain, then summed in the spatial domain over the rotation angles, to produce the perturbative correction to the current estimate of the object function. By contrast, in iODT-ISC, the perturbative Rytov phases are summed in the spatial domain over illumination angles, and then normalized in the frequency domain by the CTF, to produce the perturbative correction. Once the current estimate and its perturbative correction are available, the new estimate of the object function can be obtained for the next iteration. The algorithm produces the reconstructed object function iteratively, until the stopping condition is satisfied. Note that the first iteration of iODT-ISC will be the same as ISC-EDOF-FBPP, if the initial guess of the object function is zero (i.e., equal to the background medium).

oe-28-26-39904-i001

The normalized root-mean-square (nRMS) errors between the true fields ${u_\textrm{t}}(x = d,y;{\theta _{\textrm{ill}}})$ and the simulated forward fields $u_{\textrm{fwd}}^{(m)}(x = d,y;{\theta _{\textrm{ill}}})$ are defined as,

$${\varepsilon _\textrm{A}}(m) = \frac{{\sqrt {{{\sum\nolimits_{y,{\theta _{\textrm{ill}}}} {({|{{u_\textrm{t}}} |- |{u_{\textrm{fwd}}^{(m)}} |} )} }^2}/{N_{\textrm{pix}}}} }}{{\textrm{range}({|{{u_\textrm{t}}} |} )}},$$
$${\varepsilon _\phi }(m) = \frac{{\sqrt {{{\sum\nolimits_{y,{\theta _{\textrm{ill}}}} {[{{\textbf{\textit{W}}}(\textrm{Arg}\,{u_\textrm{t}} - \textrm{Arg}\,u_{\textrm{fwd}}^{(m)})} ]} }^2}/{N_{\textrm{pix}}}} }}{{\textrm{range}[{{\textbf{\textit{W}}}(\textrm{Arg}\,{u_\textrm{t}})} ]}},$$
where ${N_{\textrm{pix}}}$ is the total number of pixels of the sinogram, Arg returns the principal value of a phase, and ${\textbf{\textit{W}}}$ is a phase-wrapping operator. The stopping condition is based on the deviation of nRMS errors for the successive Q iterations or the maximum number of iterations,
$$\begin{array}{l} |{{\varepsilon_\textrm{A}}(m) - {\varepsilon_\textrm{A}}(m - q)} |< {\varepsilon _{\textrm{tol}}}\,\,\,\&\,\,\,|{{\varepsilon_\phi }(m) - {\varepsilon_\phi }(m - q)} |< {\varepsilon _{\textrm{tol}}},\\ \textrm{or }m > {m_{\max }}, \end{array}$$
where $q = 1,2, \ldots ,Q$, and $\varepsilon {_{\textrm{tol}}}$ is a pre-defined error tolerance.

A requirement for iODT-ISC is that the transverse dimension of the simulation window for tilted beam propagation should increase with the illumination angle, leading to an increase in the simulation time in comparison with ORC. The transverse dimension in ORC is independent of the rotation angle, whereas the dimension in ISC needs to be extended to capture the fields from tilted propagation. The simulation time for titled propagation will thus be scaled proportionally. The remaining modules in iODT-ISC are approximately as efficient as those in iODT for ORC. If the object has rotational symmetry, the data processing time can be reduced.

Although iODT-ISC, shown in Algorithm 1, is formulated in the 2D case, it can also work for three-dimensional (3D) reconstruction. The tilted beam propagation, CTF normalization, and phase unwrapping are all feasible in 3D. However, the computational time and memory requirement for 3D will increase significantly. We may improve the efficiency by use of advanced computational techniques, such as parallel programing and GPU acceleration, in the future. iODT-ISC for 2D reconstruction can be still useful, for example, for index profiling of an optical fiber based on ISC measurements, which can be more stable and faster in data acquisition. Therefore, in this paper we will focus on the 2D reconstruction by use of iODT-ISC.

3. Result

To validate the efficacy of the proposed iODT-ISC, a numerical phantom is used as the object to be imaged. The reconstruction quality is quantified by the signal-to-noise ratio (SNR),

$$SNR = 10{\log _{10}}\frac{{||{{n_\textrm{t}} - {n_\textrm{b}}} ||_{{l_2}}^2}}{{||{{n_\textrm{t}} - \widehat n} ||_{{l_2}}^2}},$$
where ${n_\textrm{t}}$ and $\widehat n$ are the RI distributions of the true object and the reconstructed object, respectively. In the examples shown in this section, the wavelength of the probing wave is ${\lambda _0} = 650$ nm, and the background RI is ${n_\textrm{b}} = 1.4584$, so that the wavelength in the background medium is ${\lambda _\textrm{b}} \approx 446$ nm. We will use this wavelength and background RI, unless otherwise specified. The phantom is sequentially illuminated by a plane wave from -60° to 60° with an increment of 1°. The true diffracted fields are computed by use of the finite difference time domain (FDTD) method [38]. In the following simulations, no constraints will be applied on the real and imaginary parts of the reconstructed object function, unless prior knowledge is specified. The field propagation solver in iODT-ISC is selected to be the tilted beam propagation method, to achieve a balance between accuracy and efficiency. L2-norm phase unwrapping [39,40] is used to remove the phase ambiguity in the imaginary part of the perturbative Rytov phase. For the stopping condition, the parameters Q and the error tolerance are chosen to be 5 and 1e-3, respectively. The simulation is implemented in MATLAB 2018a using a computer equipped with an Intel Core i7-5820K 3.3 GHz CPU and a 32GB RAM.

3.1 Two-core phantom

The object is chosen to verify the efficacy of the proposed iODT-ISC for an object with large OPD. As shown in Fig. 3(a), a two-core phantom with diameter of $40{\lambda _\textrm{b}}$ and separation of $50{\lambda _\textrm{b}}$ is created. The RI difference is $\Delta n = 0.04$. The largest OPD induced by the phantom is about 2.2λ0. The grid size in the reconstruction is $\Delta x = \Delta y = {\lambda _\textrm{b}}/5$.

 figure: Fig. 3.

Fig. 3. (a) RI distribution of a two-core phantom. Reconstructions using (b) DFI, (c) ISC-EDOF-FBPP and (d) iODT-ISC (m=18), with SNRs labelled.

Download Full Size | PDF

Figures 3(b)–3(d) show the reconstructed RI distributions using DFI, ISC-EDOF-FBPP, and iODT-ISC (m=18), respectively. Compared with the phantom RI, the SNRs of the reconstructed refractive indices are also computed and labelled. The DFI reconstruction, which is based on the focused ISC-sinogram without EDOF, has severe asymmetric features and ringing artifacts beyond the plane of best focus, as shown in Fig. 3(b). Unlike DFI, the ISC-EDOF-FBPP method or iODT-ISC (m=1) backpropagates the sinogram through the background medium, thus the depth of focus of the sinogram is extended through the object, especially for illumination angles at which two disks scatter the probing wave more independently. This is why the reconstruction in Fig. 3(c) outperforms that in Fig. 3(b), in terms of the shape. However, when the sinogram is backpropagated beyond the two-core region, the focused fields become defocused again for the remaining depth, which results in more artifacts at the edge of the reconstruction area. Furthermore, for the illumination angles at which two disks are aligned with the propagation direction (e.g., zero degree), the defocused sinogram cannot be correctly modelled by EDOF because of multiple scattering. Phase unwrapping of the complex Rytov phase is then prone to failure, which results in an underestimate of the object function. Both reasons explain why the SNR of Fig. 3(c) is compromised. To model multiple scattering, iODT-ISC forward-backward propagates the incident and measured sinograms respectively through the current estimate of the object function, instead of the background medium. The discrepancies between the forward and backward sinograms produce a more accurate perturbative Rytov phase where the 2π phase ambiguity is suppressed by phase unwrapping, and then reconstruct the perturbative object function more accurately. Therefore, compared with DFI and ISC-EDOF-FBPP, iODT-ISC (m=18) can provide the most accurate reconstruction in terms of shape, value, and artifact suppression.

The nRMS errors in the amplitude and phase of sinograms as functions of the iteration number m are plotted in Figs. 4(a) and 4(b), respectively. The nRMS error in amplitude is reduced from 7.3% to 1.8%, and the error in phase is decreased from 7.8% to 2.6%. Figures 4(c) and 4(d) show a comparison of the complex fields diffracted by the phantom (blue) and by the reconstructed object (red) after 18 iterations for the illumination angle of -40°. The size of simulation dimension is 3601×675 pixels, and the data processing time is about 10 minutes per iteration. Moreover, we have tested the two-core phantom at various contrast levels, and found that the maximum contrast iODT-ISC can adequately reconstruct is about 3.4% (Δn=0.05, OPD=2.75λ0). If the contrast continues to increase, the reconstruction will have artifacts between the two cores. However, the fields diffracted by the reconstructed object and the phantom match well, which implies that the reconstruction is trapped in a local minimum. Further increase in the contrast may result in phase-unwrapping errors and severe artifacts in the reconstruction.

 figure: Fig. 4.

Fig. 4. The nRMS errors in the amplitude (a) and phase (b) of sinograms as functions of the iteration number m. The amplitude (c) and phase (d) of the field diffracted by the phantom (blue) and by the reconstructed object (red) after 18 iterations for the illumination angle of -40°.

Download Full Size | PDF

3.2 Shepp-Logan phantom

The Shepp-Logan phantom is used to verify the efficacy of the proposed iODT-ISC in terms of resolving complex features. Here, the wavelength of the probing wave is 406 nm, and the background RI is 1.333, so that the wavelength in the background medium is about 305 nm. The phantom shown in Fig. 5(a) has a contrast of 9.3% and OPD of 0.8λ0. The grid size in the reconstruction is $\Delta x = \Delta y = {\lambda _\textrm{b}}/20$. The reconstructed RI distributions using DFI, ISC-EDOF-FBPP, iODT-ISC (m=22) are shown in Figs. 5(b)–5(d). Based on the SNR, iODT-ISC provides the best reconstruction in terms of the shape, value, and artifact suppression.

 figure: Fig. 5.

Fig. 5. (a) RI distribution of the Shepp-Logan phantom. Reconstructions using (a) DFI, (b) ISC-EDOF-FBPP, (c) iODT-ISC (m=22), with SNRs labelled.

Download Full Size | PDF

The nRMS errors in the amplitude and phase of the sinograms as functions of the iteration number m are plotted in Figs. 6(a) and 6(b). The nRMS error in amplitude is decreased from 9.4% to 1.5%, and the error in phase is reduced from 14.3% to 0.7%. The comparisons of complex fields diffracted by the phantom (blue) and by the reconstructed object (red) after 22 iterations for the illumination angle of -40° are shown in Figs. 6(c) and 6(d). The size of simulation dimension is 2561×481 pixels, and the data processing time is about 5.9 minutes per iteration.

 figure: Fig. 6.

Fig. 6. The nRMS errors in the amplitude (a) and phase (b) of sinograms as functions of the iteration number m. The amplitude (c) and phase (d) of the field diffracted by the phantom (blue) and by the reconstructed object (red) after 22 iterations for the illumination angle of -40°.

Download Full Size | PDF

3.3 Seven-core phantom

The third example is chosen to verify the efficacy of the proposed iODT-ISC for an object with large OPD and complex structures. As shown in Fig. 7(a), we employ a seven-core phantom with a contrast of 4% and OPD of 2.15λ0. The grid size in the reconstruction is $\Delta x = \Delta y = {\lambda _\textrm{b}}/4.$ The reconstructed RI distributions using DFI, ISC-EDOF-FBPP, iODT-ISC (m=39) are shown in Figs. 7(b)–7(d). The SNRs indicate that iODT-ISC provides the best reconstruction. This is also confirmed by inspection of the shape, value, and artifact suppression.

 figure: Fig. 7.

Fig. 7. (a) RI distribution of the Shepp-Logan phantom. Reconstructions using (b) DFI, (c) ISC-EDOF-FBPP, (d) iODT-ISC (m=39), with SNRs labelled.

Download Full Size | PDF

The nRMS errors in the amplitude and phase of sinograms as functions of the iteration number m are plotted in Figs. 8(a) and 8(b), respectively. The nRMS error in amplitude is lowered from 7.8% to 5.5%, and the error in phase is reduced from 11.1% to 4.7%. Figures 8(c) and 8(d) show comparisons of the complex fields diffracted by the phantom (blue) and by the reconstructed object (red) after 39 iterations for the illumination angle of -40°. The size of simulation dimension is 2873×539 pixels, and the data processing time is about 4.5 minutes per iteration.

 figure: Fig. 8.

Fig. 8. The nRMS errors in the amplitude (a) and phase (b) of sinograms as functions of the iteration number m. The amplitude (c) and phase (d) of the field diffracted by the phantom (blue) and by the reconstructed object (red) after 39 iterations for the illumination angle of -40°.

Download Full Size | PDF

4. Discussion: Missing-angle problem and constraint

As known in DFI theory, overlapping Ewald spheres are assembled in the frequency domain to form the spatial spectrum of the reconstructed object function. In ISC, the spectrum is poorly represented along the axial direction, leading to the missing-angle problem. For example, in the case of the seven-core phantom, the spectrum (base-10 log scale) shown in Fig. 9(a), is missing areas along the axial direction. This results in poor axial resolution as shown in Fig. 7(b). Prior knowledge, such as non-negativity constraint, can be iteratively employed to alleviate this missing-angle problem. If the non-negativity constraint is applied to the reconstructed object function, then the spectrum of the constrained object function can fill the missing part of the spectrum shown in Fig. 9(a). Iteratively applying the constraint in spatial domain and filling missing information in the frequency domain 30 times will form a more isotropic spectrum, as shown in Fig. 9(b). When the resultant spectrum is inverse Fourier transformed, it produces the constrained DFI result plotted in Fig. 9(c), which has higher axial resolution and increased SNR of 6.2 dB, and is superior to the unconstrained one shown in Fig. 7(b).

 figure: Fig. 9.

Fig. 9. (a) Spatial spectrum of unconstrained DFI result shown in Fig. 7(b), (b) Spatial spectrum with the missing-angle problem alleviated by use of non-negativity constraint (30 iterations), (c) refractive index distribution of constrained DFI result, with SNR labelled.

Download Full Size | PDF

The iODT-ISC results, shown in Fig. 7(d) and Fig. 8, do not rely on the non-negativity constraint to the real part of the reconstruction. We have applied the non-negativity constraint to the real part in each iteration, to further test how this constraint will help improve the reconstruction. With the non-negativity constraint enabled, iODT-ISC reconstruction converges at the 33rd iteration, which is 6 iterations faster than the unconstrained one. The SNR of the reconstruction shown in Fig. 10(a) is increased by 3.9 dB, compared with Fig. 7(d). The shape and value of each core become more accurate, while the artifact in the background region is further suppressed. As shown in Figs. 10(b) and 10(c), the nRMS error in amplitude decreases from 7.8% to 4.9%, and the error in phase reduces from 11.1% to 5.4%. FDTD is used to compute the true fields, while BPM is chosen as the propagation solver in the reconstruction for the balance of accuracy and efficiency. Therefore, the nRMS error floors are determined not only by the reconstruction, but the solver discrepancy. If the solvers are matched, the error floors in the constrained reconstruction will be lower than those in the unconstrained case. It is numerically demonstrated that prior knowledge, such as the non-negativity constraint, can improve the SNR of the reconstruction and accelerate iterative reconstruction.

 figure: Fig. 10.

Fig. 10. (a) Reconstruction using iODT-ISC (m=33) with non-negativity constraint applied; SNR is labelled. The nRMS errors in the amplitude (b) and phase (c) of sinograms as functions of the iteration number m.

Download Full Size | PDF

In DFI, the missing-angle problem is alleviated by applying a constraint to the reconstructed object function and filling missing components in the frequency domain. However, iODT-ISC with the same constraint also resolves the missing-angle problem, but in a different mechanism that compares the forward-backward propagated fields through the current reconstructed object. With the same constraint, the SNRs of the reconstructions using DFI and iODT-ISC are increased by 2.1 dB and 3.9 dB, respectively, which indicates that iODT-ISC is also better than DFI in resolving the missing-angle problem for multiply-scattering objects. Also, we have tested reconstructions from limited measurements by use of DFI and iODT-ISC for a weakly-scattering object, taking the two-core phantom as an example. All the simulation parameters are kept the same, except that the RI difference is decreased from 0.04 to 0.004. With the same non-negativity constraint, the SNRs of the reconstruction using DFI and iODT-ISC are increased by 3.5 dB and 6.3 dB, respectively. This comparison indicates that iODT-ISC also outperforms DFI in resolving the missing-angle problem for a weakly-scattering object.

5. Conclusion

As a label-free, non-invasive, and quantitative imaging technique, ODT has been used to reconstruct the refractive-index distribution of a sample from multi-view measurements in ORC or ISC. Due to its stability and fast data acquisition, ISC-based ODT has been widely used in imaging biological samples, such as cells and tissues. Conventional ODT relies on the weakly-scattering assumption, and typically fails to reconstruct multiply-scattering objects with high index contrasts, large OPDs, or complicated structures. To extend ODT to the multiply-scattering regime, iODT was developed based on iteratively matching the forward and backward propagation of the incident and diffracted fields, respectively, through the scattering area. However, iODT was originally developed for ORC and it cannot work directly for ISC-based measurements.

In this paper, we presented iODT-ISC, an update to iODT that resolved the configuration mismatch while maintaining the advantages of iODT. Compared with iODT, the update in iODT-ISC involves two aspects: tilted forward and backward beam propagation, and inversion from field discrepancies to perturbative corrections of the object function, along with the coherent transfer function normalization. Numerical simulations demonstrated that iODT-ISC can work directly for ISC and provide accurate reconstructions of samples with high contrasts, large OPDs, and complicated structures. For example, a seven-core phantom with a contrast of 4% and OPD of 2.15λ0 was accurately reconstructed. Further increase in the contrast and OPD of the object is possible with an increase in the range of illumination angles. Moreover, we have shown that under the same constraint (e.g., non-negativity), iODT-ISC, benefiting from the forward-backward field matching, can outperform the standard ODT in resolving the missing-angle problem, in terms of convergence and reconstruction quality.

Funding

National Science Foundation (1509294).

Disclosures

The authors declare no conflicts of interest.

References

1. 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(1), 266–277 (2009). [CrossRef]  

2. 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(01), 1–12 (2013). [CrossRef]  

3. M. Kujawińska, W. Krauze, A. Kus, J. Kostencka, T. Kozacki, B. Kemper, and M. Dudek, “Problems and Solutions in 3-D Analysis of Phase Biological Objects by Optical Diffraction Tomography,” Int. J. Optomechatronics 8(4), 357–372 (2014). [CrossRef]  

4. K. Kim, J. Yoon, S. Shin, S. Lee, S.-A. Yang, and Y. Park, “Optical diffraction tomography techniques for the study of cell pathophysiology,” J. Biomed. Photonics and Eng. 2(2), 020201 (2016). [CrossRef]  

5. P. Müller, M. Schürmann, and J. Guck, “ODTbrain: a Python library for full-view, dense diffraction tomography,” BMC Bioinformatics 16(1), 367–369 (2015). [CrossRef]  

6. T. Kim, R. Zhou, L. L. Goddard, and G. Popescu, “Solving inverse scattering problems in biological samples by quantitative phase imaging,” Laser Photonics Rev. 10(1), 13–39 (2016). [CrossRef]  

7. C. Hu and G. Popescu, “Quantitative Phase Imaging: Principles and Applications,” in Label-free Super-resolution Microscopy (Springer, 2019).

8. A. Kuś, W. Krauze, P. L. Makowski, and M. Kujawińska, “Holographic tomography: hardware and software solutions for 3D quantitative biomedical imaging (Invited paper),” ETRI J. 41(1), 61–72 (2019). [CrossRef]  

9. 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(6), 517–522 (2015). [CrossRef]  

10. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72(1), 156–160 (1982). [CrossRef]  

11. S. K. Debnath and Y. Park, “Real-time quantitative phase imaging with a spatial phase-shifting algorithm,” Opt. Lett. 36(23), 4677–4679 (2011). [CrossRef]  

12. W. Gorski and W. Osten, “Tomographic imaging of photonic crystal fibers,” Opt. Lett. 32(14), 1977–1979 (2007). [CrossRef]  

13. A. D. Yablon, “Multi-Wavelength Optical Fiber Refractive Index Profiling by Spatially Resolved Fourier Transform Spectroscopy,” J. Lightwave Technol. 28(4), 360–364 (2010). [CrossRef]  

14. J. Kostencka and T. Kozacki, “Computational and experimental study on accuracy of off-axis reconstructions in optical diffraction tomography,” Opt. Eng. 54(2), 024107 (2015). [CrossRef]  

15. M. Ziemczonok, A. Kuś, P. Wasylczyk, and M. Kujawińska, “3D-printed biological cell phantom for testing 3D quantitative phase imaging systems,” Sci. Rep. 9(1), 18872 (2019). [CrossRef]  

16. 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(2), 178–180 (2006). [CrossRef]  

17. M. H. B. Gilboa, Y. Roichman, and N. T. Shaked, “Tomographic phase microscopy with 180 degrees rotation of live cells in suspension by holographic optical tweezers,” Opt. Lett. 40(8), 1881–1884 (2015). [CrossRef]  

18. M. Habaza, M. Kirschbaum, C. Guernth-Marschner, G. Dardikman, I. Barnea, R. Korenstein, C. Duschl, and N. T. Shaked, “Rapid 3D Refractive-Index Imaging of Live Cells in Suspension without Labeling Using Dielectrophoretic Cell Rotation,” Adv. Sci. 4(2), 1600205 (2017). [CrossRef]  

19. F. Merola, P. Memmolo, L. Miccio, R. Savoia, M. Mugnano, A. Fontana, G. D’Ippolito, A. Sardo, A. Iolascon, A. Gambale, and P. Ferraro, “Tomographic flow cytometry by digital holography,” Light: Sci. Appl. 6(4), e16241 (2017). [CrossRef]  

20. T. Cacace, P. Memmolo, M. M. Villone, M. De Corato, M. Mugnano, M. Paturzo, P. Ferraro, and P. L. Maffettone, “Assembling and rotating erythrocyte aggregates by acoustofluidic pressure enabling full phase-contrast tomography,” Lab Chip 19(18), 3123–3132 (2019). [CrossRef]  

21. S. Shin, K. Kim, J. Yoon, and Y. Park, “Active illumination using a digital micromirror device for quantitative phase imaging,” Opt. Lett. 40(22), 5407–5410 (2015). [CrossRef]  

22. A. B. Ayoub, J. Lim, E. E. Antoine, and D. Psaltis, “Optical Diffraction Tomography Based on a Spatial Light Modulator for Biological Imaging,” in Biophotonics Congress: Optics in the Life Sciences Congress 2019 (BODA,BRAIN,NTM,OMA,OMP), The Optical Society (OSA, 2019), NS1B.6.

23. 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(13), 16933–16948 (2015). [CrossRef]  

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

25. 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(1), 59–70 (2016). [CrossRef]  

26. S. Chowdhury, M. Chen, R. Eckert, D. Ren, F. Wu, N. Repina, and L. Waller, “High-resolution 3D refractive index microscopy of multiple-scattering samples from intensity images,” Optica 6(9), 1211–1219 (2019). [CrossRef]  

27. X. Ma, W. Xiao, and F. Pan, “Optical tomographic reconstruction based on multi-slice wave propagation method,” Opt. Express 25(19), 22595–22607 (2017). [CrossRef]  

28. D. Suski, J. Winnik, and T. Kozacki, “Fast multiple-scattering holographic tomography based on the wave propagation method,” Appl. Opt. 59(5), 1397–1403 (2020). [CrossRef]  

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

30. H. 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 4(1), 73–86 (2018). [CrossRef]  

31. S. Fan, S. Smith-Dryden, J. Zhao, S. Gausmann, A. Schülzgen, G. Li, and B. E. A. Saleh, “Optical Fiber Refractive Index Profiling by Iterative Optical Diffraction Tomography,” J. Lightwave Technol. 36(24), 5754–5763 (2018). [CrossRef]  

32. S. Fan, S. Smith-Dryden, G. Li, and B. Saleh, “Reconstructing complex refractive-index of multiply-scattering media by use of iterative optical diffraction tomography,” Opt. Express 28(5), 6846–6858 (2020). [CrossRef]  

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

34. J. Kostencka, T. Kozacki, A. Kuś, B. Kemper, and M. Kujawińska, “Holographic tomography with scanning of illumination: space-domain reconstruction for spatially invariant accuracy,” Biomed. Opt. Express 7(10), 4086–4101 (2016). [CrossRef]  

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

36. A. Ritter, “Modified shifted angular spectrum method for numerical propagation at reduced spatial sampling rates,” Opt. Express 22(21), 26265–26276 (2014). [CrossRef]  

37. H.-H. Son and K. Oh, “Light propagation analysis using a translated plane angular spectrum method with the oblique plane wave incidence,” J. Opt. Soc. Am. A 32(5), 949–954 (2015). [CrossRef]  

38. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, 2005).

39. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software (Wiley, 1998).

40. J. Pineda, J. Bacca, J. Meza, L. A. Romero, H. Arguello, and A. G. Marrugo, “SPUD: simultaneous phase unwrapping and denoising algorithm for phase imaging,” Appl. Opt. 59(13), D81–D88 (2020). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (10)

Fig. 1.
Fig. 1. Configurations for ODT (a) object rotating configuration (ORC) and (b) illumination scanning configuration (ISC).
Fig. 2.
Fig. 2. Illustration of direct Fourier interpolation (DFI) for ISC.
Fig. 3.
Fig. 3. (a) RI distribution of a two-core phantom. Reconstructions using (b) DFI, (c) ISC-EDOF-FBPP and (d) iODT-ISC (m=18), with SNRs labelled.
Fig. 4.
Fig. 4. The nRMS errors in the amplitude (a) and phase (b) of sinograms as functions of the iteration number m. The amplitude (c) and phase (d) of the field diffracted by the phantom (blue) and by the reconstructed object (red) after 18 iterations for the illumination angle of -40°.
Fig. 5.
Fig. 5. (a) RI distribution of the Shepp-Logan phantom. Reconstructions using (a) DFI, (b) ISC-EDOF-FBPP, (c) iODT-ISC (m=22), with SNRs labelled.
Fig. 6.
Fig. 6. The nRMS errors in the amplitude (a) and phase (b) of sinograms as functions of the iteration number m. The amplitude (c) and phase (d) of the field diffracted by the phantom (blue) and by the reconstructed object (red) after 22 iterations for the illumination angle of -40°.
Fig. 7.
Fig. 7. (a) RI distribution of the Shepp-Logan phantom. Reconstructions using (b) DFI, (c) ISC-EDOF-FBPP, (d) iODT-ISC (m=39), with SNRs labelled.
Fig. 8.
Fig. 8. The nRMS errors in the amplitude (a) and phase (b) of sinograms as functions of the iteration number m. The amplitude (c) and phase (d) of the field diffracted by the phantom (blue) and by the reconstructed object (red) after 39 iterations for the illumination angle of -40°.
Fig. 9.
Fig. 9. (a) Spatial spectrum of unconstrained DFI result shown in Fig. 7(b), (b) Spatial spectrum with the missing-angle problem alleviated by use of non-negativity constraint (30 iterations), (c) refractive index distribution of constrained DFI result, with SNR labelled.
Fig. 10.
Fig. 10. (a) Reconstruction using iODT-ISC (m=33) with non-negativity constraint applied; SNR is labelled. The nRMS errors in the amplitude (b) and phase (c) of sinograms as functions of the iteration number m.

Equations (8)

Equations on this page are rendered with MathJax. Learn more.

F ^ ( k k inc ; θ ill ) = j 4 π k x \textit{F} ( u bgd ( y ; θ ill ) ln u ( y ; θ ill ) u bgd ( y ; θ ill ) ) ,
ϕ R ( r ; θ ill ) = ln u ( r ; θ ill ) u bgd ( r ; θ ill ) ,
f ^ ( r ) = 1 j 4 π k b θ ill ϕ R ( r ; θ ill ) Δ θ ill .
f ^ ( r ) = \textit{F} 1 [ \textit{F} ( f ^ ( r ) ) / CTF( k ) ] .
ε A ( m ) = y , θ ill ( | u t | | u fwd ( m ) | ) 2 / N pix range ( | u t | ) ,
ε ϕ ( m ) = y , θ ill [ \textit{W} ( Arg u t Arg u fwd ( m ) ) ] 2 / N pix range [ \textit{W} ( Arg u t ) ] ,
| ε A ( m ) ε A ( m q ) | < ε tol & | ε ϕ ( m ) ε ϕ ( m q ) | < ε tol , or  m > m max ,
S N R = 10 log 10 | | n t n b | | l 2 2 | | n t n ^ | | l 2 2 ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.