Open Access
Add to BibSonomy
Abstract
Optical Diffraction Tomography (ODT) is a label-free method to quantitatively estimate the 3D refractive index (RI) distributions of microscopic samples. Recently, significant efforts were directed towards methods to model multiple-scattering objects. The fidelity of reconstructions rely on accurately modelling light-matter interactions, but the efficient simulation of light propagation through high-RI structures over a large range of illumination angles is still challenging. Here we present a solution dealing with these problems, proposing a method that allows one to efficiently model the tomographic image formation for strongly scattering objects illuminated over a wide range of angles. Instead of propagating tilted plane waves we apply rotations on the illuminated object and optical field and formulate a new and robust multi-slice model suitable for high-RI contrast structures. We test reconstructions made by our approach against simulations and experiments, using rigorous solutions to Maxwell’s equations as ground truth. We find the proposed method to produce reconstructions of higher fidelity compared to conventional multi-slice methods, especially for the challenging case of strongly scattering samples where conventional reconstruction methods fail.
Published by Optica Publishing Group under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
1. Introduction
Label-free microscopy methods utilize intrinsic properties of the sample rather than fluorescent markers as contrast agents and have received considerable attention in recent years. An important example among label-free methods is quantitative phase imaging (QPI) [1], where the 2D phase map of largely transparent objects is measured with interferometry. For thin samples, the accumulated phase can be assumed to correspond approximately to the product of the refractive index (RI) and the thickness of a sample. Quantitative information about the RI is tightly linked to cell biophysical parameters such as dry mass or protein concentration [2]. In a tomographic approach, such "phase projections" are taken at different angles and fused, which yields a 3D map of the object [3]. However, for thicker samples the assumption that the phase is proportional to the integration of RI along the illumination direction no longer holds. This limits the resolution and results in artifacts in the reconstruction, such as axial elongation, as is well known [4,5].
A more sophisticated method that includes diffraction to some extent is optical diffraction tomography (ODT) [6–10]. This method enables the reconstruction of transparent micrometer-sized objects from a set of field measurements (acquiring amplitude and phase of the light scattered by the object). Multiple views of the object are obtained from which the RI distribution can be estimated. The different angular views can be obtained by scanning the illumination angle [4,11–13], rotating the object [14,15], or both [16,17]. Classical ODT is based on the Fourier Diffraction Theorem, which relates the Fourier transform of the measured scattered field to the Fourier transform of the scattering potential.
Unfortunately, classical ODT relies on the (first) Born or Rytov approximation and is therefore restricted to weakly scattering objects that cause only single-scattering events. Recently, iterative gradient-based optimization approaches have emerged to recover the RI distribution of multiple-scattering samples by modelling the nonlinear relationship between electric field and RI distribution [5,18–20]. Numerous methods have been proposed to model the non-linear relationship between field and RI or scattering potential. They can be conceptually separated into methods that solve the Lippmann-Schwinger equation and multi-slice models. The former iterate over volumes and are very accurate but are burdened by a high computational cost. The latter segment the 3D RI distribution into 2D slices and propagate light from slice to slice. These models neglect or simplify the coupling between forward and backward propagating fields, but are computationally cheap.
These methods based on numerical optimization have recently been extended to utilizing intensity-only information [21–24]. This allows for a simpler experimental setup with more cost effective illumination options such as LEDs. However, disregarding phase information may add ambiguity to the reconstruction process, which makes the optimization less stable. Therefore, a high fidelity forward model is of crucial importance to obtain a good reconstruction quality.
1.1 Current limitations
To date a stringent limitation on the above outlined ODT reconstruction from illumination under various angles is that methods to model the light-matter interaction either show a significantly reduced accuracy with increasing oblique propagation directions [5,22], high RI contrast structures [23,25] or are comparatively computationally expensive [26–29]. High tilt angles are necessary to ameliorate the missing cone problem and thus are directly connected to the quality of the reconstruction. The modeling of light propagation through high RI contrast structures is challenging, since computationally efficient methods are limited to low to moderate RI contrast structures. For a high-quality quantitative reconstruction of samples with high RI contrast one must therefore resort to more rigorous methods that are much slower.
1.2 Contributions
In the present work we present a novel approach to ODT such that both, large scan angles and high RI contrast structures can be modeled efficiently using multi-slice models. We use the following two new key ingredients for tomographic imaging:
- • We introduce a new Fast Fourier Transform (FFT)-based multi-slice method to propagate light fields through continuous high RI contrast structures. Our model represents a hybrid formulation of two popular propagation methods, the beam propagation method (BPM) [5,33] and the wave propagation method (WPM) [34,35], termed HyPM.
Moreover, the approach does not rely on the Born or Rytov approximation, which makes it a suitable method under multiple-scattering conditions.
We evaluate both modifications by numerical simulations and experiments with silica microspheres. Our results confirm a significantly higher accuracy compared to existing multi-slice methods harnessing oblique propagation. As ground truth for the simulations we use data from the infinite cylinder scattering theory [36] and full solutions of the Helmholtz equation [37], while reconstructions from experimental measurements and simulations are compared to results obtained from the Lorenz-Mie theory [3].
2. Reconstruction method
In this Section we outline the tomographic reconstruction procedure. The retrieval of the 3D RI-contrast distribution $\Delta n(\mathbf {r})$ with $\mathbf {r} = (x,y,z)^\mathrm {T}$ from a set of measurements $M_v$ under multiple-scattering represents a non-linear optimization problem. The solution $\Delta n^*$ for the RI-contrast can be written as
In Eq. (1) $S_v (\Delta n)$ represents the forward model, e.g. the simulated images, whereas $M_v$ denotes the measured images of view $v$. $V$ denotes the number of illumination angles and $||\cdot ||$ the euclidean norm. The intensity tomographic imaging problem of multiple-scattering samples represents an ill-posed non-linear inverse problem. Additionally, ODT in illumination angle scanning configuration has to deal with missing data due to the limited NA of illumination and detection optics. To ameliorate these effects we supplement Eq. (1) by a regularization term $\mathcal {R}(\Delta n) = \lambda _\mathrm {TV}\mathrm {TV}_\mathrm {I}(\Delta n)$, $\mathrm {TV}_\mathrm {I}(\Delta n)$ denoting the isotropic total variation (TV). The isotropic TV regularization can be understood as a sparsity enforcing prior and has been applied successfully in tomographic imaging. The parameter $\lambda _\mathrm {TV}$ balances the relative strength between data fidelity (first term) and regularization. To solve Eq. (1) we use proximal gradient descent [5,38–40], which is an iterative gradient-based optimization method that allows to minimize errors metrics with non-differentiable contributions.
2.1 Forward model
Our approach to model the tomographic image formation relies on the equivalence of oblique propagation and straight propagation through the rotated object followed by a rotation of the field. In Fig. 1 this correspondence is shown schematically on a 2D object. An extension to the three-dimensional case is straightforward. In (a) the light propagation is performed obliquely, whereas in (b) the propagation is carried out on the rotated object, followed by a rotation of the field. Both (a) and (b) show the same physical situation, as (b) is simply the rotated version of (a).
Fig. 1. Schematic light propagation (BPM) through an exemplary inhomogeneous RI distribution (green and blue show contour lines of the RI distribution) and a recording plane (black dashed). (a) shows the oblique propagation, (b) depicts our approach where object and field are rotated. Both (a) and (b) show the same physical setting.
2.1.1 Object rotation
Before the numerical propagation, we rotate the object represented by the RI-contrast $\Delta n (\mathbf {r})$ by a rotation matrix $\hat {\mathrm {R}}$, where the rotation is parameterized by a rotation angle $\vartheta$ and rotation axis $\mathbf {a} = (a_x,a_y,0)$, which is connected to the direction of the illumination, with $\vartheta = 0$ corresponding to normal incidence. To obtain a numerically rotated RI distribution $\Delta n_\mathrm {rot}(\mathbf {r}) = \Delta n(\hat {\mathrm {R}}^{-1} \mathbf {r})$ we use linear interpolation.
2.1.2 Propagation through high RI contrast structures
To propagate through the inhomogeneous medium we utilize a novel combination of two existing multi-slice propagation methods. One is the Beam Propagation Method (BPM) [33], which is a popular and efficient way to model the propagation of light through microscopic structures. Multi-slice models segment the RI distribution into layers positioned at planes $z = z_l$ ($l=1,\ldots,L$), where the propagation is performed sequentially between layers. It is assumed that the RI contrast within such a layer does not vary along the z dimension. Let $\mathcal {F}$, and $\mathcal {F}^{-1}$ denote the 2D Fourier and inverse Fourier transforms acting on the transverse coordinates $\mathbf {r}_\perp = (x,y)$ in position and $\mathbf {k}_\perp = (k_x,k_y)$ in spatial frequency domain. The BPM then has the following form to describe the field evolution from $E(\mathbf {r}_\perp,z_l)$ to $E(\mathbf {r}_\perp,z_{l+1})$ between the planes $z = z_l$ and $z = z_{l+1}$ through a slab of RI contrast $\Delta n(\mathbf {r}_\perp,z_l)$ of thickness $\Delta z_l = z_{l+1} - z_l$
The second method is the Wave Propagation Method (WPM) [34], which treats the the light-matter interaction more rigorously. The WPM has the following form to describe the field evolution:
The WPM does not assume small refractive index variations or paraxiality, but rather solves the one-way Helmholtz equation. It does so by considering the phase shifts arising from the propagation through the inhomogeneity separately for all positions and spatial frequency components. Therefore, the WPM described by Eq. (6) does not separate into steps made in position and spatial frequency space, like the BPM and other split-step multi-slice methods. This property allows the WPM to model high-RI contrast structures, e.g., to evaluate the performance of micro-optical components, where the propagation of light in microscopic structures up to $\Delta n = 0.5$ is correctly described [43,44].
Although being more rigorous, a shortcoming of the classical WPM is the high computational cost of evaluating Eq. (3) due to the mixing of position and spatial frequency space variables $\mathbf {r}_\perp$ and $\mathbf {k}_\perp$. However, in [43] the authors showed that Eq. (3) can be efficiently evaluated for a piece-wise constant refractive index distribution. Given a piece-wise constant structure with refractive index values $\Delta s^m$ in disjoint regions $\mathcal {J}^m$ ($m=1,\ldots,M$) and $\mathrm {K}_{l}^m(\mathbf {k}_\perp ) = \exp (i \Delta z_l k_z^m(\mathbf {k}_\perp ))$ with wave numbers $k_z^m(\mathbf {k}_\perp ) = \sqrt {(n_0+\Delta s^m)^2 k_0^2 - \mathbf {k}_\perp ^2}$, the WPM (Eq. (3)) can be evaluated by
To evaluate Eq. (3) efficiently for a continuous RI distribution, we extend Eq. (4) by a phase modulation step, incorporating the relative RI contrast with respect to the RI steps as shown in Fig. 2. The RI contrast can then be written as
The propagation step through a single slab then becomes
Eq. (6) is a combination of the BPM (Eq. (2)) and the WPM (Eq. (4)), that is a hybrid propagation method (HyPM). If only one RI step is used in Eq. (6) it reduces to the BPM, whereas for a piece-wise constant structure, it reduces to Eq. (4). Figure 2 illustrates how an exemplary distribution $\Delta n$ (black line) in (a) can be modelled by a sum of discrete RI contrast steps $\Delta s^m$ and "wrapped" distributions $\Delta n^m$ for the combined propagation method. The RI contrast steps $\Delta s^m$ can be specifically tuned for a given refractive index structure. Furthermore, Eq. (6) can be differentiated easily with respect to $\Delta n$, which is relevant for an efficient implementation in optimization frameworks. See Supplement 1, Section 1 for a detailed description of the forward propagation and and corresponding gradient operations.
Fig. 2. Combination of WPM and BPM to HyPM. (a) depicts a 1D RI contrast profile in black and discrete levels $\Delta s^m$, (b) shows the spatial regions $\mathcal {J}^m$ containing their respective "wrapped" profiles $\Delta n^m$.
The parameter $M$ describes the number of supporting points used to segment the RI-distribution. We found that it is sufficient to choose $M$ such that the supporting points $\Delta s^m$ differ by $0.04-0.05$. Furthermore, a denser sampling, e.g. a higher value for $M$, does not bring much, since the continuous phase modulation handles this range quite well. In terms of speed, every point of support $\Delta s^m$ roughly adds the time of a single inverse FFT at each propagation step. The evaluation time for object and field rotation is independent of $M$.
2.1.3 Field rotation
To match the geometry in illumination angle scanning ODT, we evaluate the field resulting from the propagation through the rotated distribution on a tilted plane, for which we utilize efficient algorithms based on full diffraction theory [31,45,46]. These approaches make use of the property of a monochromatic wave-field, where two tilted planes are related by a rotation of the respective spectra on the Ewald-sphere. Given a two dimensional scattered field $E^\mathrm {s}(\mathbf {r}_\perp,z_I)$ on a plane $\mathbf {r}=(\mathbf {r}_\perp, z_I)^\mathrm {T}$, the field at the rotated plane $\mathbf {r^\prime } = \hat {\mathrm {R}}^{-1} \mathbf {r}$ is
2.1.4 Efficiency consideration
Rotating object and field, apart from being more accurate than oblique propagation, can be much more efficient, depending on the geometry of the sample. In a concrete numerical implementation of a forward model with multi-slice propagation as described above, be it obliquely or on-axis, it is very convenient to use three-dimensional grids in the shape of cuboids. If the sample of interest approximates roughly a sphere, a computational grid with small lateral $\mathbf {r}_\perp$ sizes can be chosen if the proposed scheme is used. For such objects, oblique propagation naturally requires an extended grid in $\mathbf {r}_\perp$ to fit the transmitted scattered field.
2.1.5 Missing frequencies during rotation
By computing the field on a tilted plane and only considering the transmitted field, certain spatial frequencies are absent in the tilted spectrum, which can be seen in Fig. 3(f), as Eq. (7) is strictly true only if reflected fields in both coordinate systems, $\mathbf {r}$ and $\mathbf {r^\prime }$, are absent. However, for the samples studied in this work we found this effect to be negligible. If desired, the frequency components can be recovered by evaluating the reflected scattered field [46].
Fig. 3. Proposed tomographic image formation model: In (a) shows a rotated object (blue), through which a plane wave with an incident angle of $0^\circ$ is propagated and gives rise to a scattered field (red). In (b) we see a 2D slice of the amplitude of the resulting scattered field and (c) depicts the Fourier transform of the field in (b), which corresponds to the 2D projection of the spectrum on the Ewald sphere. (f) depicts the projection of the spectrum on the rotated Ewald sphere, which results in the field evaluated at a tilted plane shown in (e), where the incident field was added. Steps (a,b,c,f,e) correspond to the proposed approach using rotations and (d) depicts the situation for an oblique propagation.
2.1.6 Evaluation time of the forward model
The forward model outlined in this section was implemented in the Python programming language. We leverage GPU acceleration (Nvidia RTX 2060) with self written functions in PyOpenCL. To evaluate the tomographic image formation for a set of $7$ views on a grid of $250 \times 250 \times 100$, the forward model using the HyPM ($M=3$) took $60.5$ ms, from which the rotation of the field took $2.1$ ms and the rotation of the object took $5$ ms (The rotation is performed at every propagation step), as shown in Table 1 . For the BPM the evaluation took $41$ ms, with equal grid sizes as for the HyPM. The rotation of the field was evaluated on a grid of twice the size ($500 \times 500$).
Table 1. Evaluation times (in ms) for the BPM and HyPM compared to rotations of object and optical field.
2.2 Experimental setup
Figure 4 schematically depicts our tomographic imaging setup. To obtain the intensity images we use a commercial inverted microscope (Nikon Ti-E) with a custom built illumination scanning unit. A red LED (Thorlabs M625L4) serves as the illumination source, which is imaged onto a $5$ µm pinhole and subsequently collimated using an aspheric lens (Thorlabs AL2550). A $2$D galvo system (Thorlabs GVS012) driven by a Red Pitaya (STEMLab 125-14) is used with a $150$ mm achromatic lens (Thorlabs AC508-150) and an oil-immersion objective lens (Olympus UPlanFl $1.25$ NA, $60$x) in a standard $4$f configuration to create an illumination with variable angle. For imaging a water-immersion objective (Nikon CFI APO $1.15$ NA, $40$x) with the standard microscope imaging path is used, at which end a CMOS camera (Matrix Vision BF3-2071G) is attached to a port of the microscope.
3. Results
3.1 2D reconstruction from simulated data
3.1.1 Collection of cylinders
To evaluate the performance of our method we utilize the infinite-cylinder scattering theory [36] to generate artificial sinograms which are then used to reconstruct a known phantom [25,49]. Figures 5(a,e,i) show three different ground truths consisting of irregular arrangements of cylinders with diameters $d_\mathrm {cyl} = 5$, $3.75$ and $2.5$ µm, and RI contrasts of $\Delta n_\mathrm {cyl} = 0.05$, $0.075$ and $0.1$ in a medium with background refractive index $n_0 = 1$. The data are generated for an angular distribution of $91$ evenly spaced angles varying from $-45^\circ$ to $45^\circ$ and a wavelength $\lambda = 470$ nm on an evenly spaced grid with a spacing of $\Delta x= \Delta y = 100$ nm. We performed optical filtering on the generated data to remove the evanescent part of the field. Figures 5(b,f,j), (c,g,k) and (d,h,l) correspond to reconstructions using the oblique BPM, and our method of rotating object and field using the BPM and HyPM, respectively. To quantify the performance of methods we compute the error $\varepsilon$ between ground truth and reconstruction with $\varepsilon = || \Delta n - \Delta n_\mathrm {GT} ||^2 / ||\Delta n_\mathrm {GT}||^2$. The errors in Fig. 5 of the reconstructed RI distributions decrease when used with the proposed approach of rotating object and field, regardless of the multi-slice model used. The proposed HyPM shows the lowest error across all reconstructions. Increasing RI contrast, the HyPM shows the lowest error dependence and the oblique BPM configuration shows the highest dependence, with reconstructions failing for the highest RI contrast in Fig. 5(j). For a fair comparison none of the RI contrast steps $\Delta s^m$ for the HyPM corresponded to the ground truth values and we set them at equidistant intervals. The number of steps $M$ were $2$, $3$ and $4$ for the reconstructions in Fig. 5(d), (h) and (l), respectively. We tuned the regularization parameter $\lambda _\mathrm {TV}$ over a range of values for every object and propagation method and we chose the result achieving the lowest error.
Fig. 5. Tomographic reconstruction on simulated data: (a,e,i) show three different collections of 7 cylinders with diameters of $d_\mathrm {cyl} =$ $5$, $3.75$ and $2.5$ µm, and RI contrasts of $\Delta n_\mathrm {cyl} = 0.05$, $0.075$ and $0.1$, which represent the ground truths. (b,f,j) show the respective reconstructions using BPM with oblique propagation. (c,g,k) and (d,h,l) depict the reconstructions using BPM and HyPM as propagation method, respectively, with the proposed scheme of rotating object and field. Errors $\varepsilon$ of the reconstructions with respect to the true RI distributions are depicted in red boxes. The images containing the reconstructions performed with the proposed methods are contained within the blue box.
3.1.2 Reprojection error
Two FFT-based multi-slice models that have recently been proposed to model the tomographic image formation with of multiple-scattering samples are the Split-Step Non-Paraxial method [25] (SSNP) and the Multi-Layer Born method [23] (MLB). The SSNP is a non-paraxial method that propagates the field as well as the axial derivative of the field. As a split-step method it separates steps in position and spatial frequency space. The MLB separates into free-space propagation and scattering steps, which both are performed by convolutions. This method applies the Born approximation to each of the layers.
In Fig. 6 we compare the five forward models BPM, HyPM, MLB, MLR and SSNP by calculating sinograms for the assumed arrangement of cylinders shown in Fig. 5. We increased the range of illumination angles from $-80^\circ$ to $80^\circ$. The MLR (Multi-Layer Rytov) simply corresponds to the same scheme as the MLB with the exception that the Rytov instead of the Born approximation is applied at each layer. We found the Rytov approximation to be more appropriate for strongly scattering samples. However, the MLR was not stable for all propagation angles. For numerical stability, since both the SSNP and MLB/MLR models rely on convolution kernels with asymptotic behaviour for $k_z \rightarrow 0$ [23,25], we restrict its support during propagation on a region $\sqrt {\mathbf {k}_\perp ^2} < k_0 \sin (\alpha )$ with $\alpha = 85^\circ$ (MLB and SSNP) and $\alpha = 70^\circ$ (MLR). A similar procedure was employed for the rotation of the field. Each of these five multi-slice models is applied twice: in combination with oblique propagation and in combination with our proposed scheme of rotating object and field. Consequently, we obtain 30 different sinograms in total (applying 5 multi-slice methods in 2 variants to three objects), whose accuracies are evaluated by calculating the errors with respect to the sinogram set obtained from the infinite cylinder theory serving as ground truth, averaged over the three distinct cylinder distributions for each model, respectively. For the MLR method we observe instabilities for high propagation angles, if the method is used in conjunction with oblique propagation. These angles are are omitted in Fig. 6(a) and (c).
Fig. 6. Error analysis of various multi-slice propagation models compared to infinite cylinder theory serving as ground truth. In (a) the mean squared errors of the amplitudes are shown for various multi-slice methods (BPM,HyPM,MLB,MLR,SSNP) used in the oblique configuration, while (b) shows the mean squared error for the multi-slice models used in conjunction with the proposed approach, using rotations of object and field. In (c) the difference between (a) and (b) is depicted. One can see the increase in error in (a) with the propagation angle, whereas the error in (b) shows only a minor dependence on the angle with no tendency to increase. The images containing the analysis using the proposed method are contained within the blue box.
In (a) we see the mean square errors $\varepsilon _\mathrm {oblique}$ (average values for the three objects) of field amplitudes obtained with the five propagation methods used in conjunction with oblique propagation. The sinogram errors are shown in dependence of the illumination angle and increase at higher angles for all five methods. In (b) we show analogous data for the combination of the five methods with our proposed scheme of rotating object and field, where the sinogram errors $\varepsilon _\mathrm {rotations}$ show only a very small angle dependence. In (c) the difference $\Delta \varepsilon = \varepsilon _\mathrm {oblique} - \varepsilon _\mathrm {rotations}$ between (a) and (b) is depicted. Due to effects such as reflections, which are not captured by the multi-slice models, the error fluctuates strongly, even with small changes in the illumination angle.
3.1.3 Custom phantom
To further assess the performance of our proposed method, we create a numerical custom phantom, shown in Fig. 7, which was generated using a sum of spherical harmonics with random coefficients for the outer hull ($\Delta n = 0.1$) and the small inhomogeneities ($\Delta n = 0.2$) with $n_0=1$.
Fig. 7. Tomographic reconstruction on simulated phantom: (a) Ground truth, (b) depicts the reconstruction using BPM with oblique propagation, (c) the reconstruction using BPM with the proposed rotation scheme and (d) shows our proposed HyPM method with the proposed rotation scheme. The images containing the reconstructions performed with the proposed methods are contained within the blue box.
We use the modified Born Series [37] to generate the ground truth sinograms with $\lambda = 400$ nm on a grid spacing of $\Delta x= \Delta z=20$ nm. The reconstruction is performed on a coarser grid with $\Delta x= \Delta z=80$ nm. The models used for reconstruction are the same as in Fig. 5. The regularization parameter $\lambda _\mathrm {TV}$ is again tuned over a range for every method independently, and the parameter achieving the best reconstruction accuracy is chosen. Similar to Fig. 5, the reconstruction using the BPM with rotation of field and object in Fig. 7 shows a lower error with respect to BPM with oblique propagation. The error decreases further if the HyPM is used with the proposed approach using rotations, even more drastically than for the cylinder example in Fig. 5.
3.2 3D reconstruction
3.2.1 Silica sphere
For 3D reconstruction we evaluate the proposed algorithm numerically and experimentally on silica microspheres suspended in water ($n_0=1.33$). Previous work [49] has shown that a even single sphere can induce multiple-scattering, especially if the refractive index contrast is high. To obtain ground truth values for the radius and refractive index contrast, Lorenz-Mie microscopy [50] is applied to experimental measurements of a silica bead, which yielded $\Delta n = 0.109$ and $r_s=2.12$ µm. In the following we will compare the reconstructions from simulated and measured data to the retrieved values from Lorenz-Mie theory.
For simulation of the sinograms the experimental conditions were assumed and the results from simulated and experimentally measured sinograms are shown in Fig. 8 and 9. In both Figs. 8 and 9, panels (a), (b), (e) and (f) show lateral/axial $2$D sections and $1$D slices through the reconstruction for BPM with oblique propagation, while (c), (d), (g) and (h) show the respective results obtained with both of our proposed methods. The blue dashed lines in (e-h) represent the ground truth values obtained by Lorenz-Mie microscopy.
Fig. 8. Tomographic reconstructions of a 4.2 ${\mu }$m silica sphere from simulated ground truth data (Mie-theory). (a),(b),(e),(f) depict $x-y$, $x-z$, $x$, and $z$ slices of reconstructions using the conventional oblique propagation method (BPM), whereas (c),(d),(g),(h) correspond to slices of reconstructions using the proposed method. The images containing the reconstructions performed with the proposed method are contained within the blue box.
Fig. 9. Tomographic reconstruction of a 4.2 ${\mu }$m silica sphere from experimental data of a silica bead. (a),(b),(e),(f) depict $x-y$, $x-z$, $x$, and $z$ slices of reconstructions using the conventional oblique propagation method, whereas (c),(d),(g),(h) correspond to slices of reconstructions using the proposed method. The images containing the reconstructions performed with the proposed method are contained within the blue box.
The reconstructions from simulated amplitude sinograms in Fig. 8 show an overestimation of the RI value when the BPM is used with oblique propagation (a,b,e,f), whereas the proposed approaches using the HyPM in conjunction with rotation of object and field (c,d,g,h) agrees better with the ground truth values. Furthermore, the reconstructed sphere shows a higher uniformity across the sphere in Fig. 8(c,d,g,h) compared to (a,b,e,f). The reconstructions from experimental measurements in Fig. 9 show a similar result, where the BPM in the oblique configuration (a,b,e,f) overestimates the RI value and shows a similar asymmetry along z. The reconstructed bead in Fig. 9(c,d,g,h) agrees better with the values retrieved by Lorenz-Mie microscopy and exhibits higher uniformity.
4. Discussion and outlook
In this Section we summarize and contrast the results found for the various approaches studied in a multitude of simulations and experiments, which were exemplarily represented in the previous Sections.
4.1 Reconstructions
The example from cylinder-scattering theory demonstrates the higher reconstruction fidelity obtained with our proposed propagator. The BPM performed well in reconstructing multiple cylinders with moderate RI contrast ($\Delta n = 0.05$) when combined with our proposed numerical rotation step. This highlights the ability of the BPM to model the propagation through multiple-scattering samples of moderate contrast well, despite relying on paraxial approximations. If, however, used for propagation of tilted plane waves, the reconstructions with the BPM shows distortions and non-flat distributions. With increasing RI contrast, but roughly constant accumulated phase, the reconstructions are generally of decreasing fidelity, with the HyPM showing least degradation.
For high RI contrast structures ($\Delta n = 0.1-0.2$) represented, e.g., by the custom phantom generated by a random summation of spherical harmonics, the proposed HyPM shows a clear improvement over the methods using the BPM. Even when used in conjunction with our rotation approach, the BPM performs significantly poorer than the HyPM.
The comparison to Mie-theory on a spherical object provides further insight into the performance of the multi-slice scalar propagation models. The reconstructions from simulated data clearly demonstrate the advantage of the proposed rotation method over the conventional oblique propagation approach. In both simulations and measurements, the sphere reconstructions using our proposed methods agrees very well with the ground truth parameters recovered by Lorenz-Mie theory.
4.2 Reprojection errors
An interesting result is the comparison of propagation errors of the different multi-slice models, compared in the oblique and the novel implementation using rotations. Across the board, all propagation methods suffer from a growing error with increasing propagation angle in the oblique case, where non-paraxial methods such as the WPM and SSNP - not surprisingly - show the least degradation at high angles. Despite being a non-paraxial method the MLB yields the highest error. However, for samples that are less scattering the performance of the MLB is on par with other non-paraxial propagation models. When these propagation methods are included into the proposed framework of rotating object and field, the error no longer shows the tendency to increase with the angle. This indicates that the error made by the multi-slice models studied in this work, when used off axis, is at least in part inherent.
4.3 Intensity and field information
When reconstructions are performed from measurements in the form of complex fields, the axial position of the imaging plane does not affect the reconstructions, because the full amplitude and phase information is available in any plane. But, when amplitude-only information is used, the position of the imaging plane relative to the object can influence the reconstruction. A manifestation of this effect can be seen in the reconstructions of a set of cylinders or the custom phantom, visible in Fig. 5(f), Fig. 7(b) and Fig. 9(b), where the bottom part of the RI distribution, which is closer to the tilted imaging plane, tends to take on a higher value. More generally, the reconstruction adopts artifacts and asymmetries, which depend on the position of the imaging plane relative to the object. These artifacts appear strongly in the oblique propagation case and can be alleviated by using a more sophisticated method to model the tomographic image formation such as the proposed methods.
Furthermore, the presented approach of rotating object and optical field may also be applied to field-based RI retrieval.
4.4 Outlook
Furthermore, composite formulations of multi-slice models other than that of the BPM and WPM, such as SSNP and WPM or MLB and WPM could result in accurate propagation models, which can be tailored to specific samples and experimental conditions. The presented method can be extended to incorporate reflections and vectorial light properties, which for the WPM and BPM [51], as well as the field rotation [52], have been formulated. The presented approach involving rotations of object and field can thereby be incorporated into polarization sensitive models for the consideration of optically anisotropic materials [53,54].
5. Conclusion
In this work we have introduced an efficient tomographic approach to accurately reconstruct the refractive index distributions of multiple-scattering samples from intensity images. Our method is able to handle oblique illumination and high RI-contrast structures without giving erroneous results, as verified by comparing to ground truth. As such, it provides an enabling tool for situations where currently available approaches, which are comparable in terms of computational cost, would lead to erroneous predictions.
Our first insight was to reformulate the tomographic image formation using rotations on illuminated object and optical field. We showed that this approach is particularly well suited for multi-slice models. This resulted in an efficient implementation, which lowers the loss in accuracy caused by oblique propagation. The second insight was to formulate the propagation through the inhomogeneous structure as a hybrid method, which is tailored to perform accurate on-axis propagation for high RI structures.
Methods for efficient numerical propagation of light through inhomogeneous structures have applications in many fields. Apart from applications in optical imaging, the presented methods can be used for metrology and design of micro-optical components such as high-RI gradient index materials [55].
Funding
Austrian Science Fund (F6806-N36, I3984-N3).
Acknowledgement
We would like to thank Gregor Thalhammer and Nicolas Barré for fruitful discussions and valuable input.
Disclosures
The authors declare no conflicts of interest.
Data Availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. Y. Park, C. Depeursinge, and G. Popescu, “Quantitative phase imaging in biomedicine,” Nat. Photonics 12(10), 578–589 (2018). [CrossRef]
2. P. Y. Liu, L. K. Chin, W. Ser, H. F. Chen, C.-M. Hsieh, C.-H. Lee, K.-B. Sung, T. C. Ayi, P. H. Yap, B. Liedberg, K. Wang, T. Bourouina, and Y. Leprince-Wang, “Cell refractive index for cell biology and disease diagnosis: past, present and future,” Lab Chip 16(4), 634–644 (2016). [CrossRef]
3. W. Choi, C. Fang-Yen, K. Badizadegan, S. Oh, N. Lue, R. R. Dasari, and M. S. Feld, “Tomographic phase microscopy,” Nat. Methods 4(9), 717–719 (2007). [CrossRef]
4. 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]
5. 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]
6. E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Commun. 1(4), 153–156 (1969). [CrossRef]
7. A. J. Devaney, “Inverse-scattering theory within the rytov approximation,” Opt. Lett. 6(8), 374–376 (1981). [CrossRef]
8. P. Müller, M. Schürmann, and J. Guck, “The theory of diffraction tomography,” (2015).
9. D. Jin, R. Zhou, Z. Yaqoob, and P. T. C. So, “Tomographic phase microscopy: principles and applications in bioimaging [invited],” J. Opt. Soc. Am. B 34(5), B64–B77 (2017). [CrossRef]
10. V. Balasubramani, A. Kuś, H.-Y. Tu, C.-J. Cheng, M. Baczewska, W. Krauze, and M. Kujawińska, “Holographic tomography: techniques and biomedical applications [invited],” Appl. Opt. 60(10), B65–B80 (2021). [CrossRef]
11. Y. Cotte, F. Toy, P. Jourdain, N. Pavillon, D. Boss, P. Magistretti, P. Marquet, and C. Depeursinge, “Marker-free phase nanoscopy,” Nat. Photonics 7(2), 113–117 (2013). [CrossRef]
12. 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]
13. P. Zdańkowski, J. Winnik, K. Patorski, P. Gocłowski, M. Ziemczonok, M. Józwik, M. Kujawińska, and M. Trusiak, “Common-path intrinsically achromatic optical diffraction tomography,” Biomed. Opt. Express 12(7), 4219–4234 (2021). [CrossRef]
14. A. Barty, K. Nugent, A. Roberts, and D. Paganin, “Quantitative phase tomography,” Opt. Commun. 175(4-6), 329–336 (2000). [CrossRef]
15. M. Habaza, B. Gilboa, Y. Roichman, and N. T. Shaked, “Tomographic phase microscopy with 180° rotation of live cells in suspension by holographic optical tweezers,” Opt. Lett. 40(8), 1881–1884 (2015). [CrossRef]
16. B. Simon, M. Debailleul, M. Houkal, C. Ecoffet, J. Bailleul, J. Lambert, A. Spangenberg, H. Liu, O. Soppera, and O. Haeberlé, “Tomographic diffractive microscopy with isotropic resolution,” Optica 4(4), 460 (2017). [CrossRef]
17. M. Lee, K. Kim, J. Oh, and Y. Park, “Isotropically resolved label-free tomographic imaging based on tomographic moulds for optical trapping,” Light: Sci. Appl. 10(1), 102 (2021). [CrossRef]
18. L. Tian and L. Waller, “3d intensity and phase imaging from light field measurements in an LED array microscope,” Optica 2(2), 104–111 (2015). [CrossRef]
19. 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]
20. 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(6), 609–612 (2016). [CrossRef]
21. T.-A. Pham, E. Soubies, A. Goy, J. Lim, F. Soulez, D. Psaltis, and M. Unser, “Versatile reconstruction framework for diffraction tomography with intensity measurements and multiple scattering,” Opt. Express 26(3), 2749–2763 (2018). [CrossRef]
22. 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]
23. M. Chen, D. Ren, H.-Y. Liu, S. Chowdhury, and L. Waller, “Multi-layer born multiple-scattering model for 3d phase microscopy,” Optica 7(5), 394–403 (2020). [CrossRef]
24. N. Barré, R. Shivaraman, L. Ackermann, S. Moser, M. Schmidt, P. Salter, M. Booth, and A. Jesacher, “Tomographic refractive index profiling of direct laser written waveguides,” Opt. Express 29(22), 35414–35425 (2021). [CrossRef]
25. J. Lim, A. B. Ayoub, E. E. Antoine, and D. Psaltis, “High-fidelity optical diffraction tomography of multiple scattering samples,” Light: Sci. Appl. 8(1), 82 (2019). [CrossRef]
26. 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]
27. 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 4(1), 73–86 (2018). [CrossRef]
28. T. Pham, E. Soubies, A. Ayoub, J. Lim, D. Psaltis, and M. Unser, “Three-dimensional optical diffraction tomography with lippmann-schwinger model,” IEEE Trans. Comput. Imaging 6, 727–738 (2020). [CrossRef]
29. M. Lee, H. Hugonnet, and Y. Park, “Inverse problem solver for multiple light scattering using modified born series,” Optica 9(2), 177–182 (2022). [CrossRef]
30. N. Delen and B. Hooker, “Free-space beam propagation between arbitrarily oriented planes based on full diffraction theory: a fast fourier transform approach,” J. Opt. Soc. Am. A 15(4), 857–867 (1998). [CrossRef]
31. K. Matsushima, H. Schimmel, and F. Wyrowski, “Fast calculation method for optical diffraction on tilted planes by use of the angular spectrum of plane waves,” J. Opt. Soc. Am. A 20(9), 1755–1762 (2003). [CrossRef]
32. S. D. Nicola, A. Finizio, G. Pierattini, P. Ferraro, and D. Alfieri, “Angular spectrum method with correction of anamorphism for numerical reconstruction of digital holograms on tilted planes,” Opt. Express 13(24), 9935–9940 (2005). [CrossRef]
33. M. D. Feit and J. A. Fleck, “Light propagation in graded-index optical fibers,” Appl. Opt. 17(24), 3990–3998 (1978). [CrossRef]
34. K.-H. Brenner and W. Singer, “Light propagation through microlenses: a new simulation method,” Appl. Opt. 32(26), 4984–4988 (1993). [CrossRef]
35. 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]
36. J. Schäfer, S.-C. Lee, and A. Kienle, “Calculation of the near fields for the scattering of electromagnetic waves by multiple infinite cylinders at perpendicular incidence,” Journal of Quantitative Spectroscopy and Radiative Transfer 113(16), 2113–2123 (2012). [CrossRef]
37. G. Osnabrugge, S. Leedumrongwatthanakun, and I. M. Vellekoop, “A convergent born series for solving the inhomogeneous helmholtz equation in arbitrarily large media,” J. Comput. Phys. 322, 113–124 (2016). [CrossRef]
38. Y. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course (Springer Publishing Company, Incorporated, 2014), 1st ed.
39. A. Beck and M. Teboulle, “Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems,” IEEE Trans. on Image Process. 18(11), 2419–2434 (2009). [CrossRef]
40. N. Parikh, “Proximal algorithms,” Foundations and Trends in Optimization 1(3), 127–239 (2014). [CrossRef]
41. Y. Bao and T. K. Gaylord, “Clarification and unification of the obliquity factor in diffraction and scattering theories: discussion,” J. Opt. Soc. Am. A 34(10), 1738–1745 (2017). [CrossRef]
42. J. V. Roey, J. van der Donk, and P. E. Lagasse, “Beam-propagation method: analysis and assessment,” J. Opt. Soc. Am. 71(7), 803–810 (1981). [CrossRef]
43. S. Schmidt, T. Tiess, S. Schröter, R. Hambach, M. Jäger, H. Bartelt, A. Tünnermann, and H. Gross, “Wave-optical modeling beyond the thin-element-approximation,” Opt. Express 24(26), 30188–30200 (2016). [CrossRef]
44. S. Schmidt, S. Thiele, A. Herkommer, A. Tünnermann, and H. Gross, “Rotationally symmetric formulation of the wave propagation method-application to the straylight analysis of diffractive lenses,” Opt. Lett. 42(8), 1612–1615 (2017). [CrossRef]
45. T. Tommasi and B. Bianco, “Frequency analysis of light diffraction between rotated planes,” Opt. Lett. 17(8), 556–558 (1992). [CrossRef]
46. L. Onural, “Exact solution for scalar diffraction between tilted and translated planes using impulse functions over a surface,” J. Opt. Soc. Am. A 28(3), 290–295 (2011). [CrossRef]
47. D. Ruijters, B. M. ter Haar Romeny, and P. Suetens, “Efficient GPU-based texture interpolation using uniform b-splines,” Journal of Graphics Tools 13(4), 61–69 (2008). [CrossRef]
48. D. Ruijters and P. Thevenaz, “GPU prefilter for accurate cubic b-spline interpolation,” The Computer Journal 55(1), 15–20 (2010). [CrossRef]
49. J. Lim, A. Goy, M. H. Shoreh, M. Unser, and D. Psaltis, “Learning tomography assessed using mie theory,” Phys. Rev. Appl. 9(3), 034027 (2018). [CrossRef]
50. S.-H. Lee, Y. Roichman, G.-R. Yi, S.-H. Kim, S.-M. Yang, A. van Blaaderen, P. van Oostrum, and D. G. Grier, “Characterizing and tracking single colloidal particles with video holographic microscopy,” Opt. Express 15(26), 18275–18282 (2007). [CrossRef]
51. M. Fertig and K.-H. Brenner, “Vector wave propagation method,” J. Opt. Soc. Am. A 27(4), 709–717 (2010). [CrossRef]
52. S. Zhang, D. Asoubar, C. Hellmann, and F. Wyrowski, “Propagation of electromagnetic fields between non-parallel planes: a fully vectorial formulation and an efficient implementation,” Appl. Opt. 55(3), 529–538 (2016). [CrossRef]
53. A. Saba, J. Lim, A. B. Ayoub, E. E. Antoine, and D. Psaltis, “Polarization-sensitive optical diffraction tomography,” Optica 8(3), 402–408 (2021). [CrossRef]
54. S. Shin, J. Eun, S. S. Lee, C. Lee, H. Hugonnet, D. K. Yoon, S.-H. Kim, J. Jeong, and Y. Park, “Tomographic measurement of dielectric tensors at optical frequency,” Nat. Mater. 21(3), 317–324 (2022). [CrossRef]
55. C. R. Ocier, C. A. Richards, D. A. Bacon-Brown, et al., “Direct laser writing of volumetric gradient index lenses and waveguides,” Light: Sci. Appl. 9(1), 196 (2020). [CrossRef]
