Open Access
Add to BibSonomy
Abstract
Monte Carlo (MC) simulation has been widely used to study imaging procedures, including Fourier-domain optical coherence tomography (FD-OCT). Despite the broadband nature of FD-OCT, the results obtained at a single wavelength are often used in previous studies. Some wavelength-relied imaging applications, such as spectroscopic OCT (S-OCT), are unlikely to be simulated in this way due to the lack of information from the entire spectrum. Here, we propose a novel simulator for full-wavelength MC simulation of FD-OCT. All wavelengths within the emission spectrum of the light source will be simulated, and the optical properties derived from Mie theory will be applied. We further combine the inverse discrete Fourier transform (IDFT) with a probability distribution-based signal pre-processing to combat the excessive noises in the OCT signal reconstruction, which is caused by the non-uniform distribution of the scattering events at different wavelengths. Proof-of-concept simulations are conducted to show the excellent performance of the proposed simulator on signal reconstruction and optical properties extraction. Compared with the conventional method, the proposed simulator is more accurate and could better preserve the wavelength-dependent features. For example, the mean square error (MSE) computed between the backscattering coefficient extracted by the proposed simulator and the ground truth is 0.11, which is far less than the value (7.67) of the conventional method. We believe this simulator could be an effective tool to study the wavelength dependency in FD-OCT imaging as well as a preferred solution for simulating spectroscopic OCT.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Optical coherence tomography (OCT) is a noninvasive volumetric imaging technique based on low-coherence interferometry [1]. Since the first time-domain OCT (TD-OCT) system was developed, numerous derivatives have been created during the last three decades. Among them, Fourier-domain OCT (FD-OCT), first proposed by Fercher et al. in 1995 [2], is regarded as the most influential variant and has gradually replaced TD-OCT in biomedical applications. Thanks to its $\mu$m-level resolution and mm-scale penetration, FD-OCT has been deployed in various clinical settings including ophthalmic imaging [3], intracoronary guidance [4], and airway inspection [5]. Notwithstanding that FD-OCT has a wide spectrum of applications, its related theory still has room for improvement. Ongoing researches attempt to further the theoretical understanding towards fundamental topics such as phase noise [6], multiple scattering [7], image formation [8], speckle formation [9], and image reconstruction [10].
Monte Carlo (MC) method, one of the most popular numerical tools used in biophotonics, has long been acknowledged as the gold standard for theoretically investigating the light-tissue interaction [11–14]. This approach, considering light as a collection of photons, can provide the statistical optical properties of the illuminated tissue by tracking the random events of the photons experienced inside. In previous research, “MCML", the most widely used MC simulation architecture designed by Wang et al. [15], sets the foundation of most existing biological MC studies. MC method has also been applied to numerically study the imaging process of OCT [16–18]. With this method, Yao et al. [16] described a numerical approach to distinguish the contribution of the photons from a target imaging layer in the medium (Class I) and the rest of the medium (Class II) in the signal of TD-OCT. Lu et al. [17] combined the Mie theory [19] and wavelet in MC simulation to model the light transport in homogeneous turbid media. Xiong et al. [18] modeled tissue as a combination of temporarily fixed scattering particles and background medium in the simulation, to illustrate the phenomenon of exponential decay in the OCT signal. Apart from simulations of the traditional OCT imaging process, the polarization-sensitive modality [20], speckle statistic [21], importance sampling method [22], and dual-axis OCT implementation [23] were also studied with the MC method.
The aforementioned analyses all focused on TD-OCT. In recent years, there have also been studies investigating the image formation of FD-OCT based on the MC method [24–26]. However, despite of the broadband nature of FD-OCT, these existing methods still carried out the simulation with wavelength-independent parameters: only photon packets with fixed simulation parameters at a single wavelength are launched to traverse the tissue and are eventually collected for interferogram synthesis. Excessive computational cost might be a good reason that hinders MC-based methods from simulating at full wavelengths in FD-OCT [25,26]. Another possible reason is that the signal reconstruction scheme inherited from Hartinger et al. [24] is not capable of full-wavelength simulation.
In this manuscript, we propose a novel three-stage FD-OCT simulator extending the conventional single wavelength-based MC simulation to full-wavelength based on Mie theory. In the forward MC process, unlike the conventional MC method using fixed parameters such as the fixed scattering coefficient and scattering phase function (SPF) across the full spectrum, our simulator employs wavelength-dependent parameters strictly derived from Mie theory to obtain full-wavelength simulation results. For subsequent signal reconstruction, a two-stage scheme combining the IDFT with a probability distribution-based pre-processing is implemented to generate more accurate full-wavelength signals. Simulations are conducted on semi-infinite and multi-layer phantoms, which are water solutions of microspheres with different radii, demonstrating the limitation of the conventional method and the capacity of the proposed simulator on full-wavelength signal reconstruction. Optical properties extracted from the interferogram show the higher accuracy of the proposed simulator and the great potential for simulating spectroscopic OCT (S-OCT) compared with the conventional method. To our best knowledge, this is the first time to implement the full-wavelength simulation of FD-OCT based on the MC method.
2. Theory
2.1 Conventional MC-based simulation of FD-OCT
In the classical MC-based simulator of FD-OCT [24], the algorithm can be divided into two stages: (1) forward part: photon propagation, (2) inverse part: signal reconstruction.
Following the algorithm in Hartinger et al. [24], the photon propagation can be summarized into several main steps: launch photon packets, set the step size, hit boundary, reflect or transmit, detect, absorb, and scatter. In this process, the current step size $s$ is set by a single wavelength-determined scattering coefficient ${\mu _s}$ and the absorption coefficient ${\mu _a}$ subject to Beer’s law as,
where $\xi$ is a random number between 0 and 1. The next important procedure is the scattering to change the trajectory of the photon, in which the scattering angle $\theta$ is sampled from the Henyey-Greenstein (HG) SPF [27,28]. The HG function describes the probability distribution of $\theta$ as,As to the signal reconstruction stage, by copying the received signal to the entire spectrum, the current strategy synthesizes the complete sample arm field via the following equation [24] given by,
By interfering the sample signal with the reference arm fields $r({k_m})$ in a balanced detection configuration, we can generate the final spectral interferogram as,
2.2 Limitation of the conventional MC-based simulator
Since the single wavelength-based parameters, such as ${\mu _s}$ and ${p_{{\rm {HG}}}}(\theta )$, are used in the conventional simulator, photons at each wavelength are assumed to obtain the same simulation results contributing to the sample arm field. However, due to the broadband nature of FD-OCT, these optical parameters, which are largely dependent on wavelength, will change across the full spectrum [29]. Thus simulations at each single wavelength could manifest the wavelength-dependent features of FD-OCT and might lead to entirely distinct results.
It is of great significance to develop a new MC-based simulation algorithm to model FD-OCT with full-wavelength information.
2.3 Proposed full-wavelength MC-based simulator of FD-OCT
In the MC forward process, the random scattering event is the most important part with a non-absorbing medium considered. Theoretically, biological tissues could be modeled as a swarm of spherical particles that are loosely packed in the background medium [30], which can be rigorously modeled by Mie theory. Due to the lack of efficient sampling schemes [17,31], the HG function was utilized to simply approximate the Mie theory-derived SPF at one wavelength [32]. A more accurate way is to sample the scattering angles $\theta ({\lambda _m})$ with the wavelength-dependent Mie SPF described as,
For full-wavelength signal reconstruction, the OPL $l({\lambda _m})$ and diffusion reflectance $R({\lambda _m})$ acquired at each wavelength should be integrated with a full wavelength-based sample arm field synthesis formula given by,
Intuitively, by simply stacking the results at different wavelengths, we can synthesize the full-wavelength sample arm field without the extended reflectance map ${R^m}(\lambda )$ as,
This method can be viewed as the case of migating Eq. (4) to the full-wavelength situation and seems more computational-effective. While this approach is mathematically sound, it suffers from some difficulties during the numerical computation, and details will be discussed in Section 6.1.
3. Algorithm
In this section, we propose a novel full-wavelength simulator based on the MC method for FD-OCT. The simulator can be divided into three stages refined from the conventional single wavelength-based algorithm: (1) forward MC engine, (2) probability distribution-based signal pre-processing, and (3) full-wavelength signal reconstruction, as shown in Fig. 1. In the following subsections, we will illustrate the algorithm of each stage in detail.
Fig. 1. Framework of the proposed full-wavelength simulator of FD-OCT. Stage I: forward MC engine for obtaining the wavelength-dependent OPLs $l(\lambda )$ and reflectance $R(\lambda )$; Stage II: the probability distribution-based pre-processing for generating the photon-detected probability distribution $\Phi (l,\lambda )$ and the extended reflectance maps ${R^m}(\lambda )$; Stage III: full-wavelength signal reconstruction procedure for rebuilding the A-line of FD-OCT with information from the entire spectrum.
3.1 Stage I: forward MC engine
In the first stage, photons from all wavelengths are launched in the proposed simulation framework, as the flowchart shown in Fig. 2. Compared with the algorithm in Hartinger et al. [24], modifications of our algorithm are marked in red in Fig. 2(b). For each wavelength ${\lambda _m}$, independent simulation is implemented to obtain the OPLs $l({\lambda _m})$ and the diffusion reflectance $R({\lambda _m})$, which will be used for the next signal pre-processing module.
Fig. 2. (a) Exemplary photon trajectory in a multi-layer tissue; (b) Proposed forward MC simulation algorithm. ${{\rm {P}}_j}({\lambda _m})$, $j = 1, \ldots,N'$: launched photon packets at the current wavelength. The procedures different from the conventional method are marked in red in the right figure.
3.2 Stage II: probability distribution-based signal pre-processing
Directly synthesizing the spectrum interferogram by using the Eq. (10) might lead to the failure of image formation due to the non-uniform distribution of the scattering events at different wavelengths. Details will be discussed in Section 6.1. Instead, we insert a pre-processing procedure on the collected diffusion reflectance at all wavelengths, $[R_0(\lambda _0),\ldots,R_m(\lambda _m),\ldots.,R_{M-1}(\lambda _{M-1})]$. Therefore, a photon-detected probability distribution $\Phi (l,\lambda )$ and the extended two-dimensional reflectance maps ${R^m}(\lambda )$ will be generated for full-wavelength A-line reconstruction. The algorithm is in Fig. 3, which consists of the following parts:
- (1) Quantify over OPL. We first quantify the original reflectance $[R_0(\lambda _0),\ldots,R_{M-1}(\lambda _{M-1})]$ over the OPL with a predefined interval $\Delta l = \frac {{2\pi }}{{M\delta k}}$, where $\delta k$ is the linear sampling interval in the $k$-domain (see Fig. 3(a)).
- (2) Merge reflectance. $R({\lambda _m})$ in each $\Delta l$ is merged to the reflectance intensity at $M$ discrete grids $[0,\Delta l,2\Delta l, \ldots,L]$, where $L = \left ( {1 - \frac {1}{M}} \right )\frac {{2\pi }}{{\delta k}}$.
After the first two steps, we actually obtain a discrete histogram of the diffusion reflectance collected in the simulation at all wavelengths concerning the OPL, denoted as a quantified diffusion reflectance map. The map is determined by two axes: the quantified OPL and the wavelength. In the usual histogram sense, this map represents the probability of a photon being collected in the ($l(\lambda _m),\lambda _m$) coordinate, which is crucial for signal reconstruction. While it is always noisy, the following fitting and normalization steps are necessary.
- (3) Fit over wavelengths. Due to the random property of the scattering events and the finite number of launched photon packets in the MC simulation, the quantified reflectance will not be smooth enough to show the statistical characteristics of the photon packets detected. We then use the fourth-order polynomial to fit the quantified reflectance at each grid to get the fitted reflectance map.
- (4) Normalize. We normalize the values over all wavelengths to get the detected probability distribution $\Phi (l,\lambda )$ of photon packets, which implies the wavelength-dependent feature of FD-OCT.
After obtaining the photon-detected probability distribution $\Phi (l,\lambda )$, we can extrapolate a reflectance line $R({\lambda _m})$ to an extended two-dimensional reflectance map ${R^m}(\lambda )$ in two steps: locating and redistribution as shown in Fig. 3(b). The procedure is given as follows.
- (5) Locate. For a value ${R_i}({\lambda _m})$ at ${l_i}({\lambda _m})$, we locate the corresponding probability distribution curve $\Phi ({l_i}({\lambda _m}),\lambda )$. In the most cases, ${l_i}({\lambda _m})$ will not locate at quantified grids due to the high precision of the MC simulation (see Fig. 3(c)). We implement a spline interpolation between the $\Phi ({i^ - }\Delta l,\lambda )$ and $\Phi ({i^ + }\Delta l,\lambda )$ to compute the correct distribution curve $\Phi ({l_i}({\lambda _m}),\lambda )$ over wavelengths, where ${i^ - }\Delta l$ and ${i^ + }\Delta l$ are the grids before and after ${l_i}({\lambda _m})$ respectively.
- (6) Redistribute. The reflectance ${R_i}({\lambda _m})$ at one wavelength will be redistributed to all wavelengths according to $\Phi ({l_i}({\lambda _m}),\lambda )$. This operation will repeat for the whole reflectance line $R({\lambda _m})$ to finally obtain the two-dimensional reflectance map ${R^m}(\lambda )$.
Fig. 3. Algorithm of the proposed probability distribution-based signal pre-processing. (a) Procedures to generate the wavelength-dependent photon packets detection probability distribution $\Phi (l,\lambda )$; (b) Extrapolate a single wavelength reflectance line to full spectrum; (c) Method for locating probability distribution at current OPL. Red and black dots: reflectance from a scatter at one OPL.
3.3 Stage III: full-wavelength signal reconstruction
In the last stage of full-wavelength signal reconstruction, we integrate all the reflectance maps acquired from the second stage by using Eq. (9). Following Eq. (5), the full-wavelength spectral interferogram is generated and linearized in the $k$-domain. Finally, classical IDFT is applied for A-line formation.
4. Method
4.1 Simulation configuration
Water solutions of polystyrene microspheres with the RI ${n_s} = 1.58$ and different particle concentrations are used to form two phantoms mimicking the real tissue for simulations: 0.21% with the radius equal to 1 $\mu$m (Phantom1), and 0.67% with the radius equal to 2 $\mu$m (Phantom2). Detailed simulation configurations on phantom models, photon propagation, optical parameters, and photon detection are summarized in Table 1. The spectral band is chosen according to a real-world swept laser (HSL-20, Santec Inc., JP). $\mu _a$ is simply assumed as 0 in the simulation for water solutions of polystyrene microsphere since it is far less than $\mu _s$ [21,26], and the broadband scattering coefficients are designed to include some values of common-used biological tissues [9]. Only photon packets exiting from detection filters centered on the light source with the acceptance radius limited to 15-$\mu$m and the acceptance angle limited to +/- $5^{\circ }$ will be collected.
Table 1. Configuration details of simulations
The whole program is run on a high-performance computing cluster (Siyuan-1 cluster, Shanghai, China) equipped with Intel Xeon ICX Platinum 8358 (2.6 GHz, 32 cores), and the total CPU run time is about 40,000 hours.
4.2 Effect of the broadband nature of FD-OCT
We first demonstrate the effect of the broadband nature of FD-OCT in the MC-based simulation. Three sets of optical parameters at 1210 nm, 1260 nm, and 1310 nm of wavelength are chosen for simulations on Phantom1 following the conventional MC procedures as described in Section 2.1. To acquire high-quality results, a total of 204.8 billion photons are launched for each simulation. The collected diffusion reflectance together with the final A-line at each wavelength is compared to prove the limitation of the conventional single wavelength-based simulator.
4.3 Phantom imaging simulation
We then conduct the simulations on two semi-infinite phantoms (Phantom1 and Phantom2) and a two-layer phantom by using the proposed simulator.
For each simulation on the semi-infinite phantom, 200 million photon packets per wavelength are launched, which adds up to the same number of photons at the single wavelength as in the first simulation. In the case of the multi-layer phantom, 1 billion photon packets per wavelength are utilized. The photon-detected probability distributions $\Phi (l,\lambda )$ are derived in each simulation (presented as a fitted reflectance map). We later reconstruct the A-line signals from reflectance maps with Eq. (9) to show the imaging performance of the proposed simulator.
Especially, we reconstruct A-line signals at full-wavelength and center wavelength (1260 nm) for comparison on the multi-layer phantom. For the latter, we launch 1024 billion photons for the simulation, so that the input energy of both is equal.
4.4 Signal reconstruction by stacking full-wavelength simulation results
A crucial signal pre-processing (stage II) is implemented in the proposed simulator. To illustrate the necessity of this stage and the corresponding sample arm field synthesis scheme (Eq. (9)), we compare the A-line signals reconstructed with Eq. (10) and the proposed model.
4.5 Validation: optical properties extraction in the single scattering model
To validate the proposed simulator, we extract two optical properties, scattering coefficient ${\mu _s}$ and wavelength-dependent backscattering coefficient ${\mu _b}(\lambda )$, from the full-wavelength interferogram blocking the contribution of multiple scattering photons. Specifically, Phantom2 is used for simulations. For full-wavelength simulation, 200 million photon packets per wavelength are launched. Relatively, the number of photon packets adds up to 204.8 billion for the single wavelength-based simulation.
To extract the scattering coefficient, an IDFT is applied on the interferogram to obtain the single scattering signal given by Beer’s law [33,34] as
where $z$ is the penetration depth and ${\mu _a} = 0$ for all phantoms in this study. Thus, by fitting the slope of the FD-OCT signal (in logarithmic form), the scattering coefficient can be extracted. Although wavelength-dependent scattering coefficients are used in our model, the final A-line will fit a composite scattering coefficient ${\bar \mu _s}$ that is within the interval of the true values.The spectrum ${\mu _b}(\lambda )$ can be analyzed by time-frequency analysis in spectroscopic OCT (S-OCT) [35]. The spectrum ${i^2}(z,\lambda )$ in S-OCT is extracted by applying a short-time Fourier transform (STFT) to the interferogram. ${\mu _b}(\lambda )$ can be obtained by fitting the maximum value of the spectrum at each wavelength. We compare it with the Mie theory-derived function, which has been proven to be in good agreement with experiments on polystyrene microsphere solutions [36]. The function is given by
5. Results
5.1 Effect of the broadband nature of FD-OCT
We first investigate the Mie SPFs of the semi-infinite phantom (see Fig. 4(a)) at three wavelengths in Fig. 4(b). It is shown that all curves agree well in the forward scattering part (small angle), but have much diversity for the backscattering angles around $\pi$ as pointed out by the black arrows. At 1310 nm, the probability at $\pi$ angle is 43 times larger than the case at 1210 nm. Since OCT mainly collects backscattered signals, it will lead to a significant mismatch in the obtained diffusion reflectance, which can be seen in Fig. 4(c). As a consequence, the reconstructed FD-OCT A-lines at all wavelengths will be different from each other in either the intensity or the attenuation (see Fig. 4(d)). If we only use the parameters at a single wavelength such as 1210 nm for the MC process, the signal obtained at 1310 nm is greatly underestimated, which testifies to the necessity of including all wavelengths within the entire spectral band during the simulation.
Fig. 4. Comparison of the single wavelength simulation results. (a) The implementation of the infinite phantom. There is no reflection at the incident interface due to the matched RI between the background layer and the phantom; (b) Comparison of the Mie scattering phase function ${p_{{\rm {Mie}}}}(\theta ;x)$ visualized as a function of scattering angle at three wavelengths; (c) Obtained OPL-resolved diffusion reflectance; (d) Normalized A-line signals. The parameter ${\mu _s}$ is calculated as 4.1 ${\rm mm}^{-1}$ at 1210 nm, 3.7 ${\rm mm}^{-1}$ at 1310 nm, and 3.9 ${\rm mm}^{-1}$ at 1210 nm. The absorption coefficient is 0 for all simulations. Black arrows in (b): ${p_{{\rm {Mie}}}}(\pi ;x)$ values.
5.2 Phantom imaging simulation
Figure 5(a) and (c) show fitted reflectance maps along the OPL and wavelength for semi-infinite Phantom1 and Phantom2 respectively. Obviously, reflectance gradually decreases when the OPL becomes larger. When viewing the wavelength side, different fluctuation patterns are presented for two phantoms, which perceive disparate wavelength-dependent features. The extracted profiles show a closer look at this mechanism. Figure 5(b) and (d) present the normalized A-line signal by using our full-wavelength reconstruction method. A sharp peak can be found at the first interaction surface (depth = 0 mm) of the phantom in each signal, which does not clearly appear in the cases of the single wavelength (see Fig. 4(d)).
Fig. 5. OPL-resolved reflectance map and normalized A-line signal for Phantom1 and Phantom2 by the proposed simulator. (a), (c) OPL-resolved reflectance map of Phantom1 and Phantom2. (b), (d) Normalized full-wavelength A-line signal of Phantom1 and Phantom2. The profiles at 1230 nm (blue), 1290 nm (red), and center wavelength 1260 nm (black) are extracted to show the variation along the OPL direction (right panels in (a) and (c)), while the profiles at 480 $\mu$m (blue) and 160 $\mu$m (red) are extracted to show the variation over wavelength (top panels in (a) and (c)).
To further show the usefulness of the proposed simulator, we design the simulation of a two-layer phantom consisting of Phantom1 and Phantom2 as shown in Fig. 6(a). The physical thicknesses of two layers are $d_1$ = 0.3 mm (Phantom1), $d_2$ = 0.5 mm (Phantom2).
Fig. 6. Results of a two-layer phantom simulated by the proposed simulator. (a) Two-layer phantom. The physical thicknesses of Phantom1 and Phantom2 are 0.3 mm and 0.5 mm respectively. RI of background medium is fixed to 1.33; (b) OPL-resolved reflectance map for the two-layer phantom. The profiles at 1230 nm (blue), 1290 nm (red), and 1260 nm (black) are extracted to show the variation along the OPL direction (right panel), while the profiles at 160 $\mu$m (red) and 950 $\mu$m (blue) are extracted to show the variation over wavelength at different layers (top panel). (c) Normalized full-wavelength A-line signals of the two-layer phantom. ${\mu _s}(\lambda )$ values at 1260 nm are 3.9 ${\rm mm}^{-1}$ / 8.3 ${\rm mm}^{-1}$ for two layers. Black arrows: the first interface; red arrows: the second interface.
The reflectance map in Fig. 6(b) shows two clear interfaces. The locations 0.80 mm and 2.13 mm agree well with the round-trip OPL of each real interface, that is ${l_1} = 2 \times 1.33 \times 0.3 = 0.798\;{\rm { mm}}$, ${l_2} = 2 \times 1.33 \times 0.8 = 2.128\;{\rm { mm}}$, where ${l_1}$ and ${l_2}$ represent the OPLs of the first and the second interface respectively (see right panel of Fig. 6(b)). Besides, two wavelength-resolved profiles at 160 $\mu$m and 950 $\mu$m show completely different morphologies corresponding to two phantoms, respectively (see top panel of Fig. 6(b)).
Sharp interface features can also be found in full-wavelength A-lines (blue line) in Fig. 6(c). The depth positions of the three interfaces are 0 mm, 0.4 mm, and 1.06 mm, which are in good agreement with the positions of the reflectance interfaces in (b), since the depth should be half of the OPL. The two layers are well represented in the figure with the differences in both intensity and attenuation. As a comparison, the result from the conventional method (dash line) also characterizes the difference between two layers but seems to give more ambiguous boundaries between two layers.
5.3 Signal reconstruction by stacking full-wavelength simulation results
Normalized A-lines of two semi-infinite phantoms and the two-layer phantom are shown in Fig. 7. Compared to signals formed by the proposed simulator (blue lines), both A-lines reconstructed from Eq. (10) (red lines) are too noisy to present the structure of phantoms.
Fig. 7. Comparison of A-lines reconstructed by the proposed method and piling up the results of all wavelengths with Eq. (10). (a) Normalized A-line signals of semi-infinite Phantom1; (b) Normalized A-line signals of semi-infinite Phantom2; (c) Normalized A-line signals of the two-layer phantom.
5.4 Optical properties extraction in the single scattering model
After squaring and the logarithm, single scattering signals after IDFT at 1260 nm (conventional method) and full-wavelength (proposed method) are shown in Fig. 8(a) and (b). We use the first-order polynomial to fit their attenuation, and the slopes represent the scattering coefficients. The fitted ${\mu _{s,1260}}$ at 1260 nm is 7.0 ${\rm mm}^{-1}$, which deviates from the actual value of 8.3 ${\rm mm}^{-1}$ used in the simulation. As for the full-wavelength case, the inferred ${\bar \mu _s}$ is 8.2 ${\rm mm}^{-1}$, which is well fitted to the actual simulation parameter interval used ([8 ${\rm mm}^{-1}$, 8.7 ${\rm mm}^{-1}$]).
Fig. 8. Optical properties extraction. (a)(b) Extract scattering coefficients by IDFT analysis based on the conventional and the proposed method respectively; (c)(d) Extract spectra by STFT analysis based on the conventional and proposed method respectively; (e) Comparison between two fitted backscattering coefficient distributions and the theoretical function. The proposed simulator gives more accurate results.
For S-OCT analysis, an STFT with 1024 sample points and 16-point step size is applied to the interferogram. A Hann window is applied with a 256-point window size to get the final spectrum. The normalized results for the conventional method and the proposed full-wavelength simulator are shown in Fig. 8(c) and (d), respectively. We then use the fourth-order polynomial to fit the maximum values of each method over wavelengths to get ${\mu _b}(\lambda )$. Then we compare them with the theoretical result (black curve) computed from Eq. (12) in Fig. 8(e). The mean square error (MSE) is measured between the fitted ones and the ground truth. The curve obtained from the proposed simulator gives MSE as 0.11, which is far less than the value (7.67) from the conventional method-based curve.
6. Discussion
6.1 Origin of the excessive noise
Excessive noise is observed in Fig. 7 with the sample arm field synthesized by Eq. (10) (red lines), which does not appear in the results of the proposed (blue lines) and the conventional simulator (Fig. 4(d)). Here we will show that this effect comes from the non-uniform distribution of the scattering events at different wavelengths. In the conventional method, the results at a single wavelength are directly copied to the entire spectrum. The signal returned from a scatter of a certain OPL is equal along the wavelength as the toy example shown in Fig. 9(a). Then the signal will be presented as a full-spectrum interference fringe.
Fig. 9. Toy example imaging on a single scatter. (a) The conventional method to generate the interferogram; (b) Directly stack detected signals of all wavelengths to generate interferogram. Red dot: reflectance from a scatter at one OPL with energy $I_0$; (c) Spectral interferograms of the full spectrum and sub-sample spectrum for the scatter at OPL equal to 3.2 mm. A Hann window was adopted to original signals; (d) Reconstructed A-line signals by IDFT.
However, due to the randomness of scattering events in MC simulation, it is difficult to detect the signals returning at all wavelengths from the same OPL, even if a large number of photons are launched. Thus, directly stacking detected signals of all wavelengths as Eq. (10) will lead to the information loss at certain wavelengths (wavenumbers) as described in Fig. 9(b), which is equivalent to a sub-sampling interference pattern. We compare the full-spectrum and the sub-sampled interference fringes after windowing in Fig. 9(c), and apply a direct IDFT to reconstruct a single-depth A-line signal in Fig. 9(d). It is apparent that instead of an ideal Kronecker delta function (black), a spread noisy signal is obtained by the sub-sampled fringe. As the energy is distributed into the noise, the intensity at the peak becomes small. When the signals corresponding to the compressed interferograms at all OPLs are superimposed, all structures will be covered, resulting in the noise signal shown in Section 5.3.
In our work, we derive the wavelength-dependent photon detection probability distribution from simulation results at all wavelengths. Then we extrapolate the simulation results of each wavelength to the full spectrum through this probability distribution. In this way, we not only fill in the missing information to resolve the signal corruption caused by excessive noise but also retain wavelength-dependent features in both reflectance and depth information.
6.2 Computational cost
A potential limitation of the proposed simulator is the excessive computational cost. We summarize the computational time for each result with the proposed simulator in Table 2. To obtain high-quality probability distribution, an enormous amount of computational time is needed. In addition to the low computational efficiency of the MC process [23,29], the number of sample points in the spectrum also severely affects the computational time in our simulator. In the future, we will explore methods to speed up the whole algorithm and improve computational efficiency, such as the importance sampling method [22,37] and parallel computing by Compute Unified Device Architecture (CUDA) [38,39].
Table 2. Configuration details of simulations
6.3 Tradeoff between the quality of the simulated signal and the number of the photon packets
Billions of photon packets were launched in our simulations to generate OCT signals of the highest quality and the smallest error as possible. This oversampling configuration brings huge computational costs as shown in Table 2. The tradeoff between the quality of the simulated signal and photon packets should be considered. Take the simulation in Fig. 6 as an example, we change the number of launched photon packets per wavelength to 1 million, 10 million, 100 million, and 1 billion (used in Fig. 6) to show the variation of reconstructed signals in Fig. 10.
Fig. 10. The effect of the number of photon packets on the signal quality of the reconstructed two-layer phantom. (a)-(d) Reconstructed signals with photon packets launched per wavelength as 1 million, 10 million, 100 million, and 1 billion. Red dash line: the interface between two layers.
When fewer photon packets are used, results with reasonable attenuation and the correct interface might be reconstructed (see Fig. 10(c)). However, simulations with too few photon packets will lead to distorted structures as shown in Fig. 10(a) and (b), due to the errors in the photon-detected probability distribution $\Phi (l,\lambda )$. For the phantom without fine structures, the number of photons in the hundred million level might be sufficient.
6.4 Impact of anisotropy and multiple scattering on the extraction of $\mu _s$
The effective $\mu _s$ inferred by FD-OCT signals is often reduced by the anisotropy of tissue and the multiple scattering [40]. In the Mie theory-based model, the anisotropy $g$ of tissue is also considered as a wavelength-dependent quantity derived from the Mie scattering phase function, which is described as,
When the wavelength-dependent ${p_{{\rm {Mie}}}}(\theta ;{x})$ is used, the effect of the anisotropy on the extracted scattering coefficient is potentially considered.
To implement the single scattering signal model, we have already eliminated the contribution of multiple scattering photons in Section 4.5. For the simulation of Fig. 8(b), if we do not block the contribution of multiple scattering photons, the extracted scattering coefficient will be reduced to 2.8 ${\rm mm}^{-1}$ due to the multiple scattering as shown in Fig. 11.
Fig. 11. Comparison of extracted scattering coefficients with multiple scattering or without multiple scattering in Phantom1. The actual simulation parameter interval used is [8 ${\rm mm}^{-1}$, 8.7 ${\rm mm}^{-1}$].
6.5 From phantoms to biological samples
While the proposed tool performs proof-of-concept simulations on phantoms, the application of the method to the simulation of biological tissues will also have some pending issues. The first one is that there may be multiple types of particles coexisting in biological tissues, making it challenging to calculate $\mu _s(\lambda )$ for full-wavelength simulation. Using an integrating sphere to measure $\mu _s(\lambda )$ of biological tissue at different wavelengths would be a good solution. Secondly, the simulation uses an ideal pencil beam light source distant from the focused Gaussian beam commonly used in experimental optical systems [26]. How to incorporate the focusing characteristics of the beam into our simulator to reconstruct more accurate signals for biological tissues will be one of our future research topics.
7. Conclusion
To address the limitation of the conventional Monte Carlo-based simulation of FD-OCT, we propose a novel full-wavelength simulator based on Mie theory. Unlike the existing methods carrying out simulations at a single wavelength, which discards the inherent broadband nature of FD-OCT, the proposed simulator incorporates all information from the entire spectrum with the optical properties strictly derived from Mie theory. Excessive noise in the full-wavelength signal caused by missing information in the interferogram is avoided through a crucial probability distribution-guided signal synthesis scheme. Simulation results on full-wavelength reconstructed signals and extracted scattering coefficients have shown that the proposed simulator is more accurate and could better preserve the wavelength-dependent features compared with the conventional method. Moreover, the proof-of-concept simulation on backscattering coefficient extraction shows the great potential of the proposed simulator for simulating spectroscopic OCT, which is impracticable with the conventional simulator. To better understand the behavior of our simulator, we will try to compare it with other theoretical methods such as the extended Huygens-Fresnel model [40] in the future. We believe that the proposed simulator could be an effective tool to study the wavelength dependency in FD-OCT imaging as well as the theory of spectroscopic OCT.
Funding
National Natural Science Foundation of China (61905141); Open Research Fund Program of the State Key Laboratory of Low-Dimensional Quantum Physics (KF202107).
Acknowledgments
The computations in this paper were run on the Siyuan-1 cluster supported by the Center for High Performance Computing at Shanghai Jiao Tong University.
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.
References
1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef]
2. A. F. Fercher, C. K. Hitzenberger, G. Kamp, and S. Y. El-Zaiat, “Measurement of intraocular distances by backscattering spectral interferometry,” Opt. Commun. 117(1-2), 43–48 (1995). [CrossRef]
3. M. Wojtkowski, V. Srinivasan, J. G. Fujimoto, T. Ko, J. S. Schuman, A. Kowalczyk, and J. S. Duker, “Three-dimensional retinal imaging with high-speed ultrahigh-resolution optical coherence tomography,” Ophthalmology 112(10), 1734–1746 (2005). [CrossRef]
4. T. Wang, T. Pfeiffer, E. Regar, W. Wieser, H. van Beusekom, C. T. Lancee, G. Springeling, I. Krabbendam, A. F. van der Steen, R. Huber, and G. van Soest, “Heartbeat oct: in vivo intravascular megahertz-optical coherence tomography,” Biomed. Opt. Express 6(12), 5021–5032 (2015). [CrossRef]
5. H. Pahlevaninezhad, M. Khorasaninejad, Y.-W. Huang, Z. Shi, L. P. Hariri, D. C. Adams, V. Ding, A. Zhu, C.-W. Qiu, F. Capasso, and M. J. Suter, “Nano-optic endoscope for high-resolution optical coherence tomography in vivo,” Nat. Photonics 12(9), 540–547 (2018). [CrossRef]
6. Y. Ling, Y. Gan, X. Yao, and C. P. Hendon, “Phase-noise analysis of swept-source optical coherence tomography systems,” Opt. Lett. 42(7), 1333–1336 (2017). [CrossRef]
7. B. Karamata, M. Laubscher, M. Leutenegger, S. Bourquin, T. Lasser, and P. Lambelet, “Multiple scattering in optical coherence tomography. i. investigation and modeling,” J. Opt. Soc. Am. A 22(7), 1369–1379 (2005). [CrossRef]
8. Y. Ling, M. Wang, X. Yao, Y. Gan, L. Schmetterer, C. Zhou, and Y. Su, “Effect of spectral leakage on the image formation of Fourier-domain optical coherence tomography,” Opt. Lett. 45(23), 6394–6397 (2020). [CrossRef]
9. R. G. Gary, J. P. Rolland, and K. J. Parker, “Speckle statistics of biological tissues in optical coherence tomography,” Biomed. Opt. Express 12(7), 4179–4191 (2021). [CrossRef]
10. K. C. Zhou, R. Qian, A.-H. Dhalla, S. Farsiu, and J. A. Izatt, “Unified k-space theory of optical coherence tomography,” Adv. Opt. Photonics 13(2), 462–514 (2021). [CrossRef]
11. B. C. Wilson and G. Adam, “A monte carlo model for the absorption and flux distributions of light in tissue,” Med. Phys. 10(6), 824–830 (1983). [CrossRef]
12. S. T. Flock, M. S. Patterson, B. C. Wilson, and D. R. Wyman, “Monte carlo modeling of light propagation in highly scattering tissues. i. model predictions and comparison with diffusion theory,” IEEE Trans. Biomed. Eng. 36(12), 1162–1168 (1989). [CrossRef]
13. C. Zhu and Q. Liu, “Review of monte carlo modeling of light transport in tissues,” J. Biomed. Opt. 18(5), 050902 (2013). [CrossRef]
14. I. Krasnikov, A. Seteikin, and B. Roth, “Advances in the simulation of light–tissue interactions in biomedical engineering,” Biomed. Eng. Lett. 9(3), 327–337 (2019). [CrossRef]
15. L. Wang, S. L. Jacques, and L. Zheng, “MCMLl—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef]
16. G. Yao and L. V. Wang, “Monte carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44(9), 2307–2320 (1999). [CrossRef]
17. Q. Lu, X. Gan, M. Gu, and Q. Luo, “Monte Carlo modeling of optical coherence tomography imaging through turbid media,” Appl. Opt. 43(8), 1628–1637 (2004). [CrossRef]
18. G. Xiong, P. Xue, J. Wu, Q. Miao, R. Wang, and L. Ji, “Particle-fixed Monte Carlo model for optical coherence tomography,” Opt. Express 13(6), 2182–2195 (2005). [CrossRef]
19. H. C. Hulst and H. C. van de Hulst, Light Scattering by Small Particles (Courier Corporation, 1981).
20. I. Meglinski, M. Kirillin, V. Kuzmin, and R. Myllylä, “Simulation of polarization-sensitive optical coherence tomography images by a monte carlo method,” Opt. Lett. 33(14), 1581–1583 (2008). [CrossRef]
21. M. Y. Kirillin, G. Farhat, E. A. Sergeeva, M. C. Kolios, and A. Vitkin, “Speckle statistics in OCT images: Monte Carlo simulations and experimental studies,” Opt. Lett. 39(12), 3472–3475 (2014). [CrossRef]
22. I. T. Lima, A. Kalra, and S. S. Sherif, “Improved importance sampling for Monte Carlo simulation of time-domain optical coherence tomography,” Biomed. Opt. Express 2(5), 1069–1081 (2011). [CrossRef]
23. Y. Zhao, K. K. Chu, E. T. Jelly, and A. Wax, “Origin of improved depth penetration in dual-axis optical coherence tomography: a Monte Carlo study,” J. Biophotonics 12(6), e201800383 (2019). [CrossRef]
24. A. E. Hartinger, A. S. Nam, I. Chico-Calero, and B. J. Vakoc, “Monte Carlo modeling of angiographic optical coherence tomography,” Biomed. Opt. Express 5(12), 4338–4349 (2014). [CrossRef]
25. S. Zhao, “Advanced Monte Carlo simulation and machine learning for frequency domain optical coherence tomography,” Ph.D. thesis, California Institute of Technology (2016).
26. Y. Wang and L. Bai, “Accurate Monte Carlo simulation of frequency-domain optical coherence tomography,” Int. J. Numer. Method Biomed. Eng. 35(4), e3177 (2019). [CrossRef]
27. L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” The Astrophys. J. 93, 70–83 (1941). [CrossRef]
28. T. Binzoni, T. S. Leung, A. H. Gandjbakhche, D. Ruefenacht, and D. Delpy, “The use of the Henyey–Greenstein phase function in Monte Carlo simulations in biomedical optics,” Phys. Med. Biol. 51(17), N313–N322 (2006). [CrossRef]
29. L. V. Wang and H.-i. Wu, Biomedical Optics: Principles and Imaging (John Wiley & Sons, 2012).
30. L. V. Wang and H.-I. Wu, “Monte Carlo modeling of photon transport in biological tissue,” Biomedical Optics: Principles and Imaging, pp. 39–50 (2007).
31. P. Naglič, F. Pernuš, B. Likar, and M. Bürmen, “Lookup table-based sampling of the phase function for monte carlo simulations of light propagation in turbid media,” Biomed. Opt. Express 8(3), 1895–1910 (2017). [CrossRef]
32. S. L. Jacques, C. Alter, and S. A. Prahl, “Angular dependence of hene laser light scattering by human dermis,” Lasers Life Sci 1, 309–333 (1987).
33. D. J. Faber, F. J. Van Der Meer, M. C. Aalders, and T. G. van Leeuwen, “Quantitative measurement of attenuation coefficients of weakly scattering media using optical coherence tomography,” Opt. Express 12(19), 4353–4365 (2004). [CrossRef]
34. J. Yi and V. Backman, “Imaging a full set of optical scattering properties of biological tissue by inverse spectroscopic optical coherence tomography,” Opt. Lett. 37(21), 4443–4445 (2012). [CrossRef]
35. U. Morgner, W. Drexler, F. Kärtner, X. Li, C. Pitris, E. Ippen, and J. Fujimoto, “Spectroscopic optical coherence tomography,” Opt. Lett. 25(2), 111–113 (2000). [CrossRef]
36. N. Bosschaart, D. J. Faber, T. G. van Leeuwen, and M. C. Aalders Jr, “Measurements of wavelength dependent scattering and backscattering coefficients by low-coherence spectroscopy,” J. Biomed. Opt. 16(3), 030503 (2011). [CrossRef]
37. I. T. Lima, A. Kalra, H. E. Hernández-Figueroa, and S. S. Sherif, “Fast calculation of multipath diffusive reflectance in optical coherence tomography,” Biomed. Opt. Express 3(4), 692–700 (2012). [CrossRef]
38. Q. Fang and S. Yan, “Mcx cloud—a modern, scalable, high-performance and in-browser Monte Carlo simulation platform with cloud computing,” J. Biomed. Opt. 27(08), 083008 (2022). [CrossRef]
39. M. Bürmen, F. Pernuš, and P. Naglič, “Mcdataset: a public reference dataset of Monte Carlo simulated quantities for multilayered and voxelated tissues computed by massively parallel pyxopto python package,” J. Biomed. Opt. 27(8), 083012 (2022). [CrossRef]
40. L. Thrane, H. T. Yura, and P. E. Andersen, “Analysis of optical coherence tomography systems based on the extended Huygens–Fresnel principle,” J. Opt. Soc. Am. A 17(3), 484–490 (2000). [CrossRef]
