Optical tomographic reconstruction based on multi-slice wave propagation method

Xichao Ma; Wen Xiao; Feng Pan

doi:10.1364/OE.25.022595

1. Introduction

Optical tomography is a powerful technique for measuring quantitative distribution of refractive index (RI) in 3D and it is gaining increasing attentions in biomedical studies [1,2] as well as industrial inspections [3]. Its application field has been rapidly expanded [4] thanks to the significant technical advances in optical setups and reconstruction algorithms [5,6]. In optical tomography, the tested object is illuminated by a monochromatic laser beam from multiple directions by either rotating the object or scanning the illumination angle [7], and a series of transmitted complex fields are recorded holographically. A reconstruction algorithm is then implemented for recovering the 3D RI distribution, which reveals various kinds of information of objects including biomedical samples and micro-optical components [8].

The reconstruction approach is crucial to the image quality and the resolution of optical tomography. The selection of reconstruction methods depends on the adopted physical model describing the relationship between the tested RI distribution and the transmitted field. Some approximations are frequently used to linearize the relationship. In early works [9,10], the ray approximation was used in which scattering effects are neglected, making it only applicable for very weak RI contrasts. Weak scattering approximations including the first-order Born and Rytov approximations [11] also produce a linear relationship but single scattering effects are considered. Therefore, the developed method, called diffraction tomography, provides improved image quality [12]. However, if the object is thick or contains complex fast-varying structures, multiple scattering occurs and the accuracy of diffraction tomography decreases. Recently, some experimental results demonstrate that nonlinear tomographic methods are capable of improving the image quality and the resolution [13–15]. They use nonlinear forward models to describe the scattering process more elaborately and precisely. Several models have been introduced into optical tomography, such as the coupled dipole method [16], the series-expanded Born method [17] and the beam propagation method (BPM) [18]. Among these, BPM is a multi-slice approach which has been successfully used for modelling the scattering process in 3D imaging [19–21], which can also be implemented in the form of artificial neural networks [14]. In particular, the BPM-based algorithm described in [21] provides a highly promising idea to combine BPM and the optimization approach for tomographic reconstruction. BPM divides the object into a series of thin slices and propagates the light field slice-by-slice. The RI is separated into the average value and the variation, i.e. $n = \bar{n} + Δ n$ . As a result, the propagation operator is also split into two parts: propagation in a homogeneous medium $\bar{n}$ and thin-element transmission through RI variation Δn. These operations allow BPM to achieve superior results than linear models. However, its accuracy is impacted due to the slowly varying envelope approximation [22] and some other approximations when splitting the propagation operator. Therefore, BPM is mostly used for paraxial propagation and small RI variations. Its accuracy decreases when modelling light propagation through fast-varying complex structures, strong RI contrasts and in large ranges of angles [23,24]. More complex forms of BPM have been proposed for accuracy improvement [25], but the computation complexity increases significantly. The wave propagation method (WPM) was developed in 1993 to overcome the paraxial limitations of BPM [26]. Like BPM, WPM divides the object into slices. Unlike BPM, the propagation operator is not split. Instead, WPM decomposes a light field into plane-wave components and implements non-paraxial plane-wave propagation in inhomogeneous media. An integral of the propagated components yields the light field in the next slice. With more precise interpretation of light propagation, WPM provides higher accuracy when modelling light propagation in large angles, strong RI contrasts or fast-varying complex structures [24]. Therefore, WPM is frequently used in simulations for light propagation through optical components [27–29]. Although having these advantages, WPM has not been introduced into tomographic reconstruction before.

In this paper, we present a novel nonlinear tomographic reconstruction method using WPM as the forward imaging model to describe the scattering process precisely. We discuss the computational implementation of WPM. To tackle the computational complexity, we propose an efficient scheme to compute the transmitted field and its derivative. We then use an iteratively optimization method, i.e. the fast iterative shrinkage-thresholding algorithm (FISTA) [30,31], to recover the quantitative RI distribution. We also discuss the influences of the parameters in the method and how to determine their values. Experimental results demonstrate that this method possesses the capability of addressing multiple scattering problems and provides high accuracy in tomographic measurements for objects with complex fast-varying structures.

The paper is organized as follows. Section 2 provides the principles and the practical computational method of WPM. Based on this model, in Section 3, we describe the tomographic reconstruction method to make the complicated computations tractable. The key point is the efficient computation of the derivative of the transmitted fields. The main steps of the optimization approach are also presented. The experimental setup is described in Section 4. In Section 5, using experimental data, we first discuss the influences of some important parameters in the method and then demonstrate the ability of this method to address multiple scattering problems.

2. Forward imaging model

In optical tomography, the tested object is illuminated from a series of directions and the scattered light fields are collected by holographic transmission measurements. Then the RI distribution is reconstructed from these multi-angle transmitted data based on a model describing the transmission process. Therefore, the forward model is crucial in tomographic reconstruction. In this section, we introduce the principles of WPM as the forward imaging model, and provide the practical computational method.

WPM was proposed by Brenner and Singer [26] to overcome the paraxial limitations of standard BPM and was further extended into a semi-vectorial formulation [23]. Similar with BPM, WPM divides the tested object into slices of thickness Δz and propagates the light field slice-by-slice. Before the propagation through each slice, the incident field, whose complex amplitude is U(x, y, z) and propagation direction is the z direction, is decomposed into plane-wave components according to

\tilde{U} (k_{x}, k_{y}, z) = F_{x, y} {U (x, y, z)} = \iint U (x, y, z) \cdot e^{- i (k_{x} x + k_{y} y)} d x d y

where

F_{x, y}

denotes the 2D Fourier transform. The propagation through a slice is then performed according to

U (x, y, z + Δ z) = \iint \tilde{U} (k_{x}, k_{y}, z) \cdot φ (x, y, z, k_{x}, k_{y}) \cdot e^{i (k_{x} x + k_{y} y)} d k_{x} d k_{y}

φ (x, y, z, k_{x}, k_{y}) = e^{i Δ z \sqrt{k_{0}^{2} n^{2} (x, y, z + \frac{Δ z}{2}) - k_{x}^{2} - k_{y}^{2}}}

where n(x, y, z) is the inhomogeneous RI distribution, k₀ is the wavenumber of the incident field in vacuum and φ(x, y, z, k_x, k_y) is the phase shift factor. Due to the dependence of φ on both the spatial and the frequency variances, Eq. (2) cannot be implemented using an inverse Fourier transform. Instead, an integral of all the propagated plane-wave components is necessary. This makes WPM slower than BPM. However, compared with complex-formed BPMs and rigorous simulation methods such as the finite-difference time-domain method, WPM requires much less computation time while providing comparable accuracy.

WPM establishes a nonlinear relationship between the tested object n(x, y, z) and the scattered light field. By implementing non-paraxial propagation of plane-wave components in each slice, WPM takes multiple scattering into account and achieves high accuracy. Note that reflected fields are not considered in WPM as in BPM.

The integral in Eq. (2) is not suitable for numerical computation. Therefore, we implement it by a summation using a matrix form. For simplicity, the 2D version (x-z plane) of WPM is presented, while the 3D version can be written similarly. Suppose the tested region is divided into T slices (as depicted in Fig. 1). The discrete form of WPM in the t-th (t = 1, 2, …, T) slice is

U_{t} = \frac{1}{P} \sum_{p = - P / 2}^{P / 2 - 1} {\tilde{U}}_{t - 1} (p) \cdot h_{t} (p)

where

{\tilde{U}}_{t - 1} (p)

is one of the plane-wave components; p is the discrete coordinate in the frequency domain along the x direction;

U_{t} \in ℂ^{P}

is a column vector, denoting the propagated field; P is the number of pixels of each slice; the column vector

h_{t} (p) \in ℂ^{P}

is the phase shifting factor of each component with its elements being

{(h_{t} (p))}_{ξ} = e^{i Δ z \sqrt{k_{0}^{2} n_{t} {(ξ)}^{2} - {(p Δ k_{x})}^{2}}} \cdot e^{i (ξ Δ x \cdot p Δ k_{x})}, ξ = - P / 2, \dots, P / 2 - 1

where n_t(ξ) is the RI distribution in the t-th slice; ξ is the discrete coordinate in the spatial domain along the x direction; Δx and Δk_x are the resolutions in the spatial and the frequency domains, respectively.

Fig. 1 Schematic diagram of WPM.

Download Full Size | PDF

Equation (4) is clumsy to implement. Therefore, we provide an alternative form using matrix multiplication, which has the potential to achieve improved efficiency when implemented in frequently-used numerical computing software such as MATLAB. Equation (4) is rewritten as

U_{t} = \frac{1}{P} H_{t} F U_{t - 1}

where column vectors

U_{t - 1}, U_{t} \in ℂ^{P}

are the light fields before and after propagation, respectively;

F \in ℂ^{P \times P}

is the discrete Fourier transform matrix so that FU_t_-1 produces the plane-wave components;

H_{t} \in ℂ^{P \times P}

is formed as

H_{t} = {[h_{t} (- \frac{P}{2}), h_{t} (- \frac{P}{2} + 1), ..., h_{t} (\frac{P}{2} - 1)]}_{P \times P}

The matrix multiplication in Eq. (6) efficiently implements the summation of all the propagated components. Therefore, from a known incident field U₀, we can iteratively obtain the propagating field in each slice as well as the exit field U_T.

We tested the WPM model in simulation as shown in Fig. 2. The domains in the x and z directions are both [-8.12 μm, 8.12 μm] with a sample interval of 145 nm. The tested object in Fig. 2(a) possesses fast-varying structures. The disks (diameter: 1.45 μm, RI: 1.46) are placed evenly in the medium (RI: 1.4) with an interval of 0.29 μm. A plane wave (wavelength: 532 nm) propagates along the z direction. The beam amplitude computed using WPM is shown in Fig. 2(b). As a reference, the Mie scattering theory [32] is applied to provide the ground truth and the computed amplitude of exit waves using different methods are compared in Fig. 2(c). This comparison indicates that WPM achieves higher accuracy in the investigation of complex-structured objects.

Fig. 2 Propagation of a plane wave in simulation. (a) The object with fast-varying complex structures. (b) Beam amplitude computed with WPM. (c) Comparison of exit wave amplitude computed with WPM, BPM and the Mie theory.

Download Full Size | PDF

3. Tomographic reconstruction

In optical tomography, the result is obtained from the recorded multi-angle data using a certain reconstruction method. Due to the nonlinear nature of WPM, the transmitted data cannot be directly mapped into the object space. Instead, an iterative approach is required. The reconstruction method used here is an optimization scheme, which interprets the task as the following minimization problem:

\hat{n} = \underset{n \in Γ}{\arg \min} Γ {Φ (n)} = \underset{n \in Γ}{\arg \min} {D (n) + τ \cdot R (n)}

where τ > 0 is the regularization parameter which controls the amount of regularization and the vector

n \in ℝ^{(P \times T)}

contains the RI distribution. The set

Γ \subseteq ℝ^{(P \times T)}

enforces some constraints on the RI values such as the non-negativity constraint. In the cost function Φ(n), the data-fidelity term D(n) and the regularization term R(n) are as follows:

D (n) = \frac{1}{N} \sum_{j = 1}^{N} D^{(j)} (n) = \frac{1}{2 N} \sum_{j = 1}^{N} {‖ U_{E}^{(j)} - U_{T}^{(j)} ‖}_{2}^{2}

R (n) = \sum_{j = 1}^{P \times T} \sqrt{{(\nabla_{x} n)}_{j}^{2} + {(\nabla_{z} n)}_{j}^{2}}

where N is the number of transmission measurements from different angles; D⁽^j⁾(n) denotes the quadratic error term of a single measurement;

U_{E}^{(j)}

and

U_{T}^{(j)}

are the experimental data and the computed transmitted field of an angle, respectively; ∇_x and ∇_z compute the discrete gradients along the x and z directions, respectively. The regularizer here is the frequently used isotropic TV which preserves the edges and reduces noise. Solving the problem in Eq. (8) generally depends on gradient descent methods. A main difficulty is that the nonlinear nature of the WPM model makes it inconvenient to compute the gradient of D(n). Therefore, we propose an efficient scheme to compute the derivative of the transmitted field, which is critical for the gradient of D(n). Then the popular FISTA is applied for minimization of the cost function.

3.1 Computation of gradient

In the cost function Φ(n), a quadratic term is introduced to evaluate the error between the experimental data and the estimated light field of each rotation angle, written as

D^{(j)} (n) = \frac{1}{2} {‖ U_{E}^{(j)} - U_{T}^{(j)} ‖}_{2}^{2}, j = 1, 2, \dots, N

where

U_{E}^{(j)} \in ℂ^{P}

is the experimentally recorded complex field and

U_{T}^{(j)} \in ℂ^{P}

is the transmitted field computed with WPM. D⁽^j⁾(n) considers a single rotation angle.

An important part of the reconstruction algorithm is computing the gradient of D(n), which is the average of the gradients of all D⁽^j⁾(n). In the following part, we discard the index (j) and describe a slice-by-slice scheme to compute the gradient for one angle. The gradients for other angles are computed in the same means.

The gradient is computed according to

\nabla D (n) = (\frac{\partial D (n)}{\partial n_{1}}, \frac{\partial D (n)}{\partial n_{2}}, \dots, \frac{\partial D (n)}{\partial n_{P \times T}}) = Re {{(U_{T} - U_{E})}^{H} \frac{\partial U_{T}}{\partial n}}

where ^H denotes the Hermitian transpose. In this equation, the most important part is the computation of the derivative of U_T with respective to n. This vector derivative gives a P × (P × T) matrix formed as

\frac{\partial U_{T}}{\partial n} = {[\frac{\partial U_{T}}{\partial n_{1}}, \frac{\partial U_{T}}{\partial n_{2}}, \dots, \frac{\partial U_{T}}{\partial n_{P \times T}}]}_{P \times (P \times T)}

where the columns are

\frac{\partial U_{T}}{\partial n_{j}} = {[\frac{\partial {(U_{T})}_{1}}{\partial n_{j}}, \frac{\partial {(U_{T})}_{2}}{\partial n_{j}}, \dots, \frac{\partial {(U_{T})}_{P}}{\partial n_{j}}]}^{T r}, j = 1, 2, \dots, P \times T

where ^Tr is the transpose operator.

Equation (4) provides an iterative approach to the derivative, written as

\frac{\partial U_{t}}{\partial n} = \frac{1}{P} \sum_{p = - P / 2}^{P / 2 - 1} [\frac{\partial h_{t} (p)}{\partial n} \cdot {\tilde{U}}_{t - 1} (p) + h_{t} (p) {(F \frac{\partial U_{t - 1}}{\partial n})}_{p th row}], t = 1, 2, \dots, T

where ∂h_t(p)/∂n is also a P × (P × T) matrix formed as

\frac{\partial h_{t} (p)}{\partial n} = {[\begin{matrix} O_{P \times P} & \dots & diag (b_{t} (p)) & \dots & O_{P \times P} \end{matrix}]}_{P \times (P \times T)}

where O denotes the zero matrices and diag(v) creates a square diagonal matrix with the elements of the vector v on the main diagonal. The vector

b_{t} (p) \in ℂ^{P}

contains the derivative of h_t(p) with the elements being

\begin{matrix} {(b_{t} (p))}_{ξ} = \frac{\partial {(h_{t} (p))}_{ξ}}{\partial n_{t} (ξ)} = \frac{i Δ z k_{0}^{2} n_{t} (ξ)}{\sqrt{k_{0}^{2} n_{t} {(ξ)}^{2} - {(p Δ k_{x})}^{2}}} e^{i (Δ z \sqrt{k_{0}^{2} n_{t} {(ξ)}^{2} - {(p Δ k_{x})}^{2}} + ξ Δ x \cdot p Δ k_{x})} \\ ξ = - P / 2, \dots, P / 2 - 1 \end{matrix}

The matrix ∂h_t(p)/∂n consists of T square matrices and Eq. (16) indicates that the t-th one is a diagonal matrix while the others are zero matrices.

Equation (15) is clumsy, making the computation time-consuming. Therefore, we provide another scheme using matrix multiplication, which can improve the computational efficiency in typical numerical computing software. Equation (15) can be rewritten as

\frac{\partial U_{t}}{\partial n} = \frac{1}{P} (J_{t} + H_{t} F \frac{\partial U_{t - 1}}{\partial n})

J_{t} = {[\begin{matrix} O_{P \times P} & \dots & diag (B_{t} F U_{t - 1}) & \dots & O_{P \times P} \end{matrix}]}_{P \times (P \times T)}

B_{t} = {[b_{t} (- \frac{P}{2}), b_{t} (- \frac{P}{2} + 1), \dots, b_{t} (\frac{P}{2} - 1)]}_{P \times P}

where

J_{t} \in ℂ^{P \times (P \times T)}

is a block matrix consisting of T square matrices of which the t-th one is diagonal and the others are zero matrices;

B_{t} \in ℂ^{P \times P}

is a square matrix with the columns corresponding with different plane-wave components.

Since the incident light field U₀ does not depend on the RI distribution n, we have

\frac{\partial U_{0}}{\partial n} = O_{P \times (P \times T)}

Therefore, in order to obtain the gradient ∇D(n), we should first compute the complex light field propagating through the current RI distribution and keep the intermediate fields in all slices. Then using Eq. (18) and the initial condition Eq. (21), we can iteratively compute the derivative of U_T and obtain ∇D(n) with Eq. (12). Note that some additional operations are necessary in the object-rotating configuration of tomography. Before computing the propagating field, the current RI distribution should be rotated by the corresponding angle while the beam propagates in a fixed direction. After obtaining the gradient, ∇D(n) should be rotated reversely by the same angle. In contrast, in the illumination-scanning configuration, the RI distribution remains still but the beams propagate in different angles, thus avoiding the rotations. The average of the gradients for all rotation angles is then used in the reconstruction algorithm to minimize the cost function Φ(n).

3.2 Iterative reconstruction

In order to minimize the cost function Φ(n), FISTA is used and the main steps are presented in Fig. 3. This framework is inspired by [21]. Due to similar schemes of slice-by-slice computations of BPM and WPM, the iterative reconstruction can be implemented in a similar manner. By computing ∇D(n), FISTA iteratively renews the RI distribution so that the cost function value decreases. The Fast Gradient Projection (FGP) algorithm in [31] is used for TV regularization with the maximum number of sub-iterations being 20. The step size ρ is 10⁻⁵ in this paper. A step size which shrinks during iteration can also be used as in [21]. The iteration stops when a maximum number of iteration k_max is reached or the relative difference between the solutions of two successive iterations is smaller than a tolerance δ (typically 10⁻⁴).

Fig. 3 Flow chart of the reconstruction algorithm.

Download Full Size | PDF

Some key parameters influence the performance of the method, such as the initial guess, the imposed constraints of RI (denoted with Γ) and the regularization parameter τ. Generally, these parameters are not determined by theoretical rules but by practical results. The influences of them are discussed in Section 5.1.

In the algorithm, the most time-consuming part is computing ∇D(n), which requires to compute ∇D⁽^j⁾(n) for all N rotation angles. However, in practice, N is large (usually hundreds), making the computation highly time-consuming. Therefore, a stochastic version is preferred, of which the original idea was described in [21]. Instead of computing N gradients, only a small number ( $N^{'} < N$ ) of gradients are computed. The rotation angles are randomly permuted and divided into groups of $N^{'}$ items (the last group may contain less than $N^{'}$ items). In the iterations, these groups are successively selected and the gradients of $N^{'}$ angles are computed, their average being ∇D(n). After the last group is used, all the angles are randomly permuted and divided again. This selection process guarantees that each angle is used with equal probability. In this way, the computation time can be significantly reduced while the results are similar. Further, the parallel computing technique can be adopted to compute multiple gradients simultaneously. $N^{'} = 5$ is used in the experiments in the remainder of the paper.

4. Experimental setup

In order to prove the validity of the proposed method, we recorded experimental data. The experimental setup is based on a Mach-Zehnder off-axis holographic interferometer in the object-rotating configuration, as depicted in Fig. 4. The laser beam (λ = 532 nm) is divided into an object beam and a reference beam, both of which are filtered and expanded by spatial filters and collimated to produce plane waves. The tested object, an optical fiber here, is held vertically on a micro-positioning stage and immersed in the RI-matching liquid (glycerol, RI: 1.454) in the cuvette. The object beam traverses the object and is captured by an imaging system composed of the MO (Sigma Koki, 20 × , NA = 0.4, infinity corrected) and the tube lens (f = 200 mm). The interference pattern is recorded by the camera (Lu125M, Lumenera Co., Ontario, Canada; 1024 × 1024; pixel pitch: 6.7 μm) as a digital hologram. 180 holograms are recorded with a 2° step during a full rotation of 360°.

Fig. 4 Schematic diagram of the experimental setup: M, mirror; PBS, polarizing beam splitter; λ/2, half-wave plate; SF, spatial filter with collimating lens; MPS, micro-positioning stage with the measured sample; C, cuvette containing RI-matching liquid; MO, microscope objective; TL, tube lens; BS, beam splitter; OB, object beam; RB, reference beam.

Download Full Size | PDF

5. Experimental results

Using experimental data, we first discuss the influences of some important parameters and determine the proper selections of them, including the initial guess, the RI constraint and the regularization parameter. A proper initial guess provides a good start of the iteration process. The constraint sets boundaries to the RI values during iteration, avoiding wrong directions of convergence. The regularization parameter controls the amount of TV regularization, thus influencing the image quality. Then we demonstrate the performance of the method in dealing with multiple scattering problems by comparing the results of this WPM-based algorithm (WPMA) with those of filtered backpropagation (FBP, a linear method) [33] and the BPM-based algorithm (BPMA) [21].

5.1 Discussion of the influences of parameters

First, we measured a four-core optical fiber (cladding diameter: 125 μm), whose structure is relatively simple. The RI difference between the cladding and the medium is about 0.003. In this condition, the weak scattering approximation can be regarded to be valid. Therefore, the reconstruction using FBP can be used as a reference. The reconstructed images contain 128 × 128 pixels with a sample interval of 1.34 μm.

The influences of different initial guesses are first discussed. The reconstructions are demonstrated in Fig. 5. The non-negativity constraint is imposed and the regularization parameter is 0.7. Two kinds of initial guesses are used: a constant-value (medium RI) initialization and a warm initialization which is the result of FBP. Figures 5(b) and 5(c) show that both initializations lead to similar results, which agree with the FBP reconstruction. In Fig. 5(d), the cost function value indicates the error between the experimental data and the computed data from the reconstruction. It is seen that the warm initialization gives much smaller initial error, but both curves converge to similar values except that the final error of FBP initialization is a little smaller. Although both initializations perform well, the reconstruction of a linear method as a warm initialization is preferred because it provides faster convergence and relatively more accurate RI values. Moreover, for complex structures, a warm initialization generally guarantees the algorithm to converge to the global solution and produce a correct result, while a constant initialization may lead to a local minimum and give a wrong result.

Fig. 5 Reconstructions of the four-core fiber with different initial guesses in WPMA. (a) Reconstruction using FBP as a reference. (b, c) Reconstructions using WPMA initialized with the FBP result and a constant value (medium RI), respectively. (d) Cost function (τ = 0.7) value (logarithmic scale) as the function of iteration number. Scale bar: 30 μm.

Download Full Size | PDF

Next, we use the same object to investigate the influences of different constraints. The results shown in Fig. 6 are initialized with a constant value (medium RI) and the regularization parameter is 0.7. The non-negativity constraint is frequently adopted in tomographic reconstruction [10,12]. It is based on the knowledge that the RI of the object is not smaller than that of the medium. If no constraint is introduced, the estimated RI is not bounded. In Fig. 6, it is clear that WPMA provides a correct result with the non-negativity constraint but a wrong result without a constraint. This indicates that without a constraint and a proper initialization, the algorithm is very likely to be trapped in a local minimum and produce an incorrect result. If a warm initialization is adopted, the constraints are not necessary (as in Fig. 7) because when the initial guess is close to the final result, the values are not likely to reach the boundaries during iteration.

Fig. 6 Reconstructions of the four-core fiber with different constraints. (a) Reconstruction using FBP as a reference. (b, c) Reconstructions using WPMA with non-negativity constraint and no constraint, respectively. (d) Cost function (τ = 0.7) value (logarithmic scale) as the function of iteration number. Scale bar: 30 μm.

Download Full Size | PDF

Fig. 7 Reconstructions of the PCF with different amounts of regularization. (a-c) Reconstructions using WPMA with τ = 0, 0.1 and 3, respectively. (d) A cross section of PCF imaged with a scanning electron microscope (SEM). Scale bar: 60 μm.

Download Full Size | PDF

The parameter τ controls the amount of TV regularization, which can reduce noise and provide smooth results. However, as a disadvantage, the best value of τ cannot be decided theoretically. It depends on the size and complexity of the object. In general, τ can be determined in the range [10⁻², 10¹]. We use another object to investigate τ, which is a photonic crystal fiber (PCF) (cladding diameter: 245 μm). All the air holes are filled with the RI-matching liquid so that the RI in them should be the same with the medium. The RI difference between the cladding and the medium is about 0.004. The reconstructed images contain 256 × 256 pixels and the sample interval is 1.34 μm. The FBP result is used as the initial guess and no constraint is imposed. The reconstructions with different amounts of regularization are shown in Fig. 7. The result with τ = 0 (Fig. 7(a), no regularization) suffers from strong noise so that the RI values fluctuate heavily. In contrast, in the result with τ = 3 (Fig. 7(c), excessive regularization), the structures deform and the air holes shrink. If an even larger τ is used, the holes will be blended together. Comparing with the SEM image in Fig. 7(d), the parameter τ = 0.1 proves to be a balance between the reduction of noise and the fidelity of the reconstruction.

5.2 Reconstruction of complex-structured objects

After the proper selections of the parameters are determined, in this section, the experimental data are reconstructed with different methods, including FBP, BPMA and WPMA. The object is also the PCF, which is a typical kind of objects that possess complex fast-varying structures, inducing severe multiple scattering effects. In BPMA and WPMA, the reconstruction of FBP is used as the initial guess and no RI constraint is imposed. The regularization parameter of WPMA is 0.1 according to Fig. 7. The same parameter of BPMA is determined to be 0.15.

The reconstructions are shown in Fig. 8. The linear FBP method results in distortions of structures, especially in the central area, as in Fig. 8(a). Both BPMA and WPMA can correct for multiple scattering so that their results in Figs. 8(b) and 8(c) surpass the FBP result. Further, the zoom-in images indicate that WPMA provides the best result. The air holes even in the innermost layer are distinct and of regular shapes while in the BPMA result, there are more deformations and the central area protrudes. In Fig. 8(d), both algorithms converge well but the curve of WPMA falls faster. These results indicate that although both BPMA and WPMA possess the capability of dealing with multiple scattering problems, WPMA provides better results and less distortions.

Fig. 8 Reconstructions of the PCF. (a-c) Results of FBP, BPMA and WPMA, respectively. (d) Cost function value as a function of iteration number. Scale bar: 60 μm.

Download Full Size | PDF

The computation times of these three methods are compared in practice. The algorithms are implemented in MATLAB R2016b and the computations are performed on a desktop computer equipped with Intel Xeon E5-2670 2.3 GHz CPU and 64 GB RAM. Averagely, for a reconstruction image of 128 × 128 pixels, FBP costs only 0.23 s. BPMA and WPMA need 11.4 s and 18.2 s, respectively, per iteration, or 38 min and 60 min, respectively, for 200 iterations. For a 256 × 256 image, FBP needs 0.47 s. BPMA costs 198 s per iteration and WPMA 307 s. These comparisons indicate that WPMA costs approximately 60% more time than BPMA does. This comes from the higher complexity of WPM, but their efficiencies are still comparable.

In order to reach a situation where the multiple scattering effects are even more severe, the RI-matching liquid in the experiment is diluted with water, leaving the RI difference between the cladding of the PCF and the medium (also the air holes) approximately 0.008. The regularization parameters of WPMA and BPMA are 0.1 and 0.2, respectively. The reconstructions are shown in Fig. 9. The result of FBP in Fig. 9(a) suffers from serious distortions, while the results of BPMA and WPMA in Figs. 9(b) and 9(c) are much better due to their ability of correcting for multiple scattering. Between these two methods, WPMA provides clearer structures and less deformations. These results demonstrate that WPMA provides better results than BPMA when dealing with severe multiple scattering problems.

Fig. 9 Reconstructions of the PCF in diluted medium. (a-c) Results of FBP, BPMA and WPMA, respectively. (d) Cost function value as the function of iteration number. Scale bar: 60 μm.

Download Full Size | PDF

6. Conclusions

In summary, this paper describes a novel nonlinear tomographic reconstruction method using WPM as the forward model. As an accurate multi-slice model, WPM has not been introduced into tomographic reconstruction before. Efficient schemes are proposed to compute both the propagating light fields and their derivatives with respective to the RI distribution. Using a popular algorithm, FISTA, the cost function is iteratively minimized so that the algorithm converges to a result that replicates the true RI distribution. Using the experimental data of complex-structured objects, we demonstrate that WPMA possesses the capability to address multiple scattering problems and provides more accurate results than BPMA under severe multiple scattering.

In practical applications, it is worth noting that the image quality and the accuracy of WPMA depends on the selections of some important parameters, including the initial guess, the constraints of RI and the regularization parameter. Using experimental data, we discuss the influences of these parameters. The initial guess provides a starting point of the iterative algorithm and the RI constraints bound the directions of convergence. If improper selections are used, the algorithm is very likely to be trapped in a local minimum thus yielding a wrong result. With a warm initialization, generally the result of a linear method, the global minimum and a correct reconstruction can be successfully reached. In this situation, the RI constraint is not necessary. The regularization parameter τ is a key parameter influencing the accuracy of WPMA. With an underestimated value, the reconstruction will suffer from noise and other disturbances. If τ is excessive, the smoothing effect will be exaggerated, erasing some detailed structures. In general, the optimal value of τ should be determined according to practical results. By choosing proper selections of the parameters, we can obtain accurate reconstructions using WPMA.

A disadvantage of this method is its computational cost. Not being able to use Fast-Fourier-Transforms, WPM is slower than BPM and more time is needed computing the gradients of light fields. However, the computation time of WPMA is not prohibitive. Repeated comparisons show that the time of WPMA is approximately 60% more than that of BPMA, making an acceptable amount of time. Besides, like in BPMA, the number of rotation angles used in each iteration can be substantially reduced in the stochastic version (from 180 to 5). For many objects whose transmission images of different directions are similar, this number can even be reduced to 1, forming an online-learning fashion. Sharing this advantage with BPMA, WPMA can provide comparable computation speeds. Using some advanced computing techniques, such as parallel computing and GPU acceleration, the computational efficiency can be further improved.

This method is of significance in tomographic measurements of complex-structured objects and thick objects inducing multiple scattering. Although only optical fibers are measured in the object-rotating configuration in this paper, this proposed method can also be used in the illumination-scanning configuration and for various kinds of objects including biomedical samples. The experimental results show that nonlinear methods are promising to further improve the image quality and the resolution of optical tomography.

References and links

1. K. Kim, S. Lee, J. Yoon, J. Heo, C. Choi, and Y. Park, “Three-dimensional label-free imaging and quantification of lipid droplets in live hepatocytes,” Sci. Rep. 6(1), 36815 (2016). [CrossRef] [PubMed]

2. H. Park, S. Lee, M. Ji, K. Kim, Y. Son, S. Jang, and Y. Park, “Measuring cell surface area and deformability of individual human red blood cells over blood storage using quantitative phase imaging,” Sci. Rep. 6(1), 34257 (2016). [CrossRef] [PubMed]

3. K. Kim, J. Yoon, and Y. Park, “Large-scale optical diffraction tomography for inspection of optical plastic lenses,” Opt. Lett. 41(5), 934–937 (2016). [CrossRef] [PubMed]

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

5. 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–463 (2017). [CrossRef]

6. 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] [PubMed]

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

8. 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 Eng. 2, 20201 (2016).

9. F. Charrière, N. Pavillon, T. Colomb, C. Depeursinge, T. J. Heger, E. A. D. Mitchell, P. Marquet, and B. Rappaz, “Living specimen tomography by digital holographic microscopy: morphometry of testate amoeba,” Opt. Express 14(16), 7005–7013 (2006). [CrossRef] [PubMed]

10. 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] [PubMed]

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

12. 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] [PubMed]

13. G. Maire, F. Drsek, J. Girard, H. Giovannini, A. Talneau, D. Konan, K. Belkebir, P. C. Chaumet, and A. Sentenac, “Experimental demonstration of quantitative imaging beyond Abbe’s limit with optical diffraction tomography,” Phys. Rev. Lett. 102(21), 213905 (2009). [CrossRef] [PubMed]

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

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

16. K. Belkebir, P. C. Chaumet, and A. Sentenac, “Influence of multiple scattering on three-dimensional imaging with optical diffraction tomography,” J. Opt. Soc. Am. A 23(3), 586–595 (2006). [CrossRef] [PubMed]

17. H.-Y. Liu, D. Liu, H. Mansour, P. Boufounos, L. Waller, and U. S. Kamilov, “SEAGLE: Robust Computational Imaging under Multiple Scattering,” in Imaging and Applied Optics 2017, OSA Technical Digest (online) (Optical Society of America, 2017), paper MM4C.1.

18. M. D. Feit and J. A. Fleck Jr., “Light propagation in graded-index optical fibers,” Appl. Opt. 17(24), 3990–3998 (1978). [CrossRef] [PubMed]

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

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

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

22. G. L. Pedrola, Beam Propagation Method for Design of Optical Waveguide Devices (Wiley, 2015), Chap. 2.

23. M. Fertig and K.-H. Brenner, “Vector wave propagation method,” J. Opt. Soc. Am. A 27(4), 709–717 (2010). [CrossRef] [PubMed]

24. 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] [PubMed]

25. E. V. Bekker, P. Sewell, T. M. Benson, and A. Vukovic, “Wide-angle alternating-direction implicit finite-difference beam propagation method,” J. Lightwave Technol. 27(14), 2595–2604 (2009). [CrossRef]

26. K.-H. Brenner and W. Singer, “Light propagation through microlenses: a new simulation method,” Appl. Opt. 32(26), 4984–4988 (1993). [CrossRef] [PubMed]

27. W. Singer, M. Testorf, and K.-H. Brenner, “Gradient-index microlenses: numerical investigation of different spherical index profiles with the wave propagation method,” Appl. Opt. 34(13), 2165–2171 (1995). [CrossRef] [PubMed]

28. J. Kostencka, T. Kozacki, M. Dudek, and M. Kujawińska, “Noise suppressed optical diffraction tomography with autofocus correction,” Opt. Express 22(5), 5731–5745 (2014). [CrossRef] [PubMed]

29. 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] [PubMed]

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

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

32. 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,” J. Quant. Spectrosc. Radiat. Transf. 113(16), 2113–2123 (2012). [CrossRef]

33. A. J. Devaney, “A filtered backpropagation algorithm for diffraction tomography,” Ultrason. Imaging 4(4), 336–350 (1982). [CrossRef] [PubMed]

Optical tomographic reconstruction based on multi-slice wave propagation method

Abstract

1. Introduction

2. Forward imaging model

3. Tomographic reconstruction

3.1 Computation of gradient

3.2 Iterative reconstruction

4. Experimental setup

5. Experimental results

5.1 Discussion of the influences of parameters

5.2 Reconstruction of complex-structured objects

6. Conclusions

References and links

Cited By

Figures (9)

Equations (21)

Optics Express