Open Access
Add to BibSonomy
Abstract
The Monte Carlo (MC) method is one of the most widely used numerical tools to model the light interaction with tissue. However, due to the low photon collection efficiency and the need to simulate the entire emission spectrum, it is computationally expensive to simulate the full-spectrum backscattered diffuse reflectance (F-BDR). Here, we propose an acceleration scheme based on importance sampling (IS). We derive the biasing sampling function tailored for simulating BDR based on the two-term scattering phase function (TT). The parameters of the TT function at different wavelengths are directly obtained by fitting the Mie scattering phase function. Subsequently, we incorporate the TT function and its corresponding biased function into the redefined IS process and realize the accelerated simulation of F-BDR. Phantom simulations based on the Fourier-domain optical coherence tomography (FD-OCT) are conducted to demonstrate the efficiency of the proposed method. Compared to the original simulator without IS, our proposed method achieves a 373× acceleration in simulating the F-BDR of the multi-layer phantom with a relative mean square error (rMSE) of less than 2%. Besides, by parallelly computing A-lines, our method enables the simulation of an entire B-scan in less than 0.4 hours. To our best knowledge, it is the first time that a volumetric OCT image of a complex phantom is simulated. We believe that the proposed acceleration method can be readily applied to fast simulations of various F-BDR-dependent applications. The source codes of this manuscript are also publicly available online.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
An accurate analysis of the wavelength-dependent backscattered light from biological tissues is highly desired, because it can potentially provide non-invasive monitoring of tissue structure [1], microscopic characteristics [2], and spectroscopic information [3]. A precise analytical or numerical modeling of full-spectrum backscattered diffuse reflectance (F-BDR) is crucial in elucidating the relationship between tissue parameters and imaging results, which could help us gain a better understanding of various phenomena including speckle [4,5], dispersion [6,7], and phase noise [8,9].
Monte Carlo (MC) method is a widely used technique for numerically simulating light-tissue interactions in biological tissues [10–12]. Combined with the accurate Mie scattering theory [13], the MC method has been applied to a variety of biophotonics studies, including polarized light transport [14,15], diffuse reflectance spectroscopy [16,17], fluorescence [18,19], and optical coherence tomography (OCT) [20,21]. Despite its effectiveness in simulating photon transport, the relatively low photon collection efficiency of backscattered photons makes it computationally expensive in the simulation of F-BDR. This issue is further exacerbated by the wavelength-dependent nature of light-tissue interactions, as simulating each wavelength independently results in additional computational costs [21].
Various techniques are proposed to address this issue [22–26]. Among them, an variance reduction technique called importance sampling (IS) [27] can be deployed with commonly used CPUs. It has been widely applied in diverse fields, including confocal microscopy [22], OCT [24,25], and optical wireless communication [28]. The essence of IS for simulating BDR lies in constructing a biased scattering function based on the normal scattering phase function (SPF) that is dominated by forward scattering. Subsequently, each biased event is weighted by the likelihood ratio of observing it in an unbiased simulation [22,24,25,28].
Unfortunately, the IS method cannot be directly applied to the simulations which uses Mie SPF in non-closed form. It is because the IS requires an analytic SPF for deriving the biased scattering function. Therefore, a closed-form SPF approximating the Mie function is needed. The classical approximation utilized in most IS-based methods is the Henyey-Greenstein (HG) SPF [29,30]. However, it has been demonstrated to be inaccurate for simulating the backscattered events [31,32]. Additionally, it does not incorporate wavelength dependency, which makes it unsuitable for simulating the F-BDR.
To address the above issues, in this manuscript we propose a novel acceleration scheme for the MC simulation of F-BDR based on IS. A two-term five-parameter SPF (TT) proposed by Jacques and McCormick [32] is utilized to approximate the Mie SPF. Specifically, we first derive the TT-biased function for the TT-based IS process of MC simulation to bias photons towards the detector. Then, by fitting the TT function and the wavelength-dependent Mie function, we obtain the wavelength-dependent parameters of TT and its biased form. Finally, the aforementioned functions are incorporated into the IS-based scattering event to achieve the accelerated simulation of F-BDR. The method is validated on MC simulations of polystyrene microsphere phantoms using Fourier-domain OCT (FD-OCT). Normalized F-BDR for the multi-layer phantom shows that the proposed IS-accelerated method can increase the speed of the original simulator by 373 times with the relative mean square error (rMSE) less than 2%. As a result, the reduced computational cost allows us to complete the B-mode simulation within 0.4 hours and makes it possible to present the FD-OCT volumetric simulation result of a complex phantom. To our best knowledge, this is the first time to achieve a fast simulation of F-BDR via MC method. The proposed IS-based acceleration scheme could act as a useful tool for studies of light-tissue interactions. The source codes are now publicly available online at [33].
2. Theory
2.1 Two-term scattering phase function
The TT function introduced by Jacques and McCormick [32] is a linear combination of two one-term SPFs derived from the generating function for Gegenbauer polynomials, which is also known as the GK function [34]. This one-term SPF can be expressed as follows,
2.2 Biased two-term SPF for simulating backscattered diffuse reflectance
In [25], Lima et al. introduced an importance sampling technique for simulating BDR. This method involves manually biasing the scattered photons towards the detector and subsequently sampling an angle $\theta _B$ by a biased SPF. The $\theta _B$ is defined as the angle between the new propagation direction and the detector direction, with a limit of $\cos \theta _B \in [0,1]$. However, this method is developed on the basis of HG function and is inherently not suitable for simulating F-BDR.
Inspired by this approach, we derive the biased form of the TT function, where we also limit the sampled biased angle to less than $\pi /2$. This ensures that after the biased scattering, the photon will not keep moving away from the detector, which improves the detection efficiency. The original cosine of the scattering angle in a GK function is described as [32],
It could be verified that this function satisfies the normalization requirement presented in Eq. (3), where the upper limit of the integral is modified to $\pi /2$. It is worth noting that when $\alpha _B=1/2$, Eq. (8) reduces to Eq. (2) in [25]. Consequently, we extend Eq. (8) to two directions with the procedure discribed in Eq. (4), thereby deriving the TT-biased SPF (TT-B), ${p_{{\rm {TT}\hbox{-}{\rm B}}}}(\theta _B ;{{\textbf {U}}_{B}} )=p_{\rm {TT}\hbox{-}{\rm B}}(\theta _B,g_{B,f},\alpha _{B,f},g_{B,b},\alpha _{B,b},C_{B})$.
2.3 Calculation of wavelength-dependent parameters in TT
Since the TT function does not inherently incorporate the wavelength dependency, we resort to Mie SPF which accurately models the wavelength-dependent scattering [13,36] as we can theoretically model biological tissues as a collection of spherical particles sparsely distributed in a background medium. The Mie SPF is designed to factor in particle parameters and the wavelength of the electromagnetic wave, which can be expressed as,
Hence, it is plausible to fit the TT function of unknown parameters and the Mie function at distinct wavelengths, thereby yielding the wavelength-dependent TT functions. These TT functions can be used to derive the corresponding TT-B function for accelerated simulations while ensuring the wavelength dependency of the simulated F-BDR.
3. Algorithm
In this section, we introduce a novel IS-based acceleration algorithm for the MC simulation of F-BDR using the TT function and its biased form TT-B function. The proposed algorithm consists of two steps, as illustrated in Fig. 1. In the first step, we use the five-parameter TT functions to fit Mie functions at distinct wavelengths, thus preserving the wavelength-dependent nature of scattering (see Fig. 1(a)). In the second step, we incorporate the TT and TT-B functions into an IS method, resulting in a biased sampling of backscattering angles, as shown in Fig. 1(b).
Fig. 1. Procedures of the proposed IS-based acceleration algorithm. (a) Step I: fit the TT function with unknown parameters and a determined Mie function at the wavelength of $\lambda$; (b) Step II: incorporate TT and TT-B functions into the biased sampling of the scattering event. The photon will be biased towards the detector; (c) A brief workflow of MC simulation with IS; (d) The algorithm of the scattering with IS.
3.1 Fitting the five-parameter TT function to the Mie function by nonlinear least square method
The nonlinear least square method (NLSM) [37] is utilized to fit the TT function and the Mie function by optimizing the parameters in the TT function, as demonstrated in Fig. 1(a). Specifically, the parameters of the TT function at wavelength $\lambda$ is ${{\textbf {U}}_\lambda } = \{ {g_{f,\lambda }},{g_{b,\lambda }},{\alpha _{f,\lambda }},{\alpha _{b,\lambda }},{C_{\lambda }}\}$, and the optimization process for each wavelength follows the formula as,
3.2 Importance sampling based on the TT function
To accelerate the MC simulation of F-BDR, we apply an IS method that biases the photons towards the detector in scattering events as shown in Fig. 1(c). This method incorporates the wavelength-dependent TT functions and TT-B functions as shown in Fig. 1(b) and Fig. 2. The flow chart for determining the scattering angle using the IS method is presented in Fig. 1(d). During the photon movement, the moving direction is defined as $\textbf{u}=( u_x, u_y,u_z)$, while the direction of the detector is $\textbf{v}=( v_x, v_y,v_z)$, and the new moving direction after scattering is $\textbf{u'}=({u}'_x,{u}'_y,{u}'_z)$. For the simulation at each wavelength $\lambda$, supposing all the photons are shot vertically from the detector into the tissue, that is $u_z=1$. The main parts of IS are listed below:
- A. First scattering event bias. When the current photon has not experienced the biased scattering and is moving away from the detector, which means ${u}_z>0$, it will experience the first biased scattering (First-bias), as shown in Fig. 2(a). For the biased angle $\theta _B=\arccos (\textbf{u'} \cdot \textbf{v})$ sampled from ${p_{{\rm {TT}\hbox{-}{\rm B}}}}$ and the corresponding “true” scattering angle $\theta _S=\arccos (\textbf{u} \cdot {\textbf{u'}})$, the likelihood ratio is calculated as,
- B. Split.We use a photon splitting mechanism similar to that in [25]. When a photon has been first biased scattered to ${\textbf{u'}}$, another unbiased photon will be generated at the original position with a large likelihood ratio $L' = 1-L({\theta _B};\lambda )$, and then be scattered by the angle $\theta$ sampled from ${p_{{\rm {TT}}}}(\theta ;{{\textbf {U}}_\lambda } )$, as shown in Fig. 2(b). We repeatedly sample this scattering angle until it is less than 90 degrees to avoid the unbiased photon being directly scattered back to the detector.
- C. Additional scattering bias. When the current photon is moving away from the detector (${u}_z>0$) after the first biased scattering, it will experience an additional biased scattering (Add-bias), as shown in Fig. 2(c). Thus the likelihood ratio is calculated as,$$L({\theta _B};\lambda) = \frac{{{p_{{\textrm{TT}}}}(\theta_S ;{{\textbf{U}}_\lambda} )}}{{a{p_{{\textrm{TT-B}}}}(\theta_B ;{{\textbf{U}}_{\lambda}} ) + (1 - a){p_{{\textrm{TT}}}}(\theta_S ;{{\textbf{U}}_\lambda} )}}$$where $a$ is the probability to control whether the biased scattering happens.
Fig. 2. Exemplary illustration of IS. (a) First-bias event when $u_z>0$; (b) Split the photon after the first-bias event; (c) When $u_z>0$, apply the add-bias after the first-bias.
After the scattering with IS, all likelihood ratios are multiplied to $L_{\rm all}$. The weight $W$ of the current photon will be modified to $W' = L_{\rm all}W$ before it is collected.
4. Method
To validate the efficacy of our proposed IS-accelerated simulation algorithm, we apply it to the simulation of FD-OCT. Specifically, we employ the proposed method to speed up the full-spectrum simulation for FD-OCT, following the strategy proposed in [21].
4.1 Simulation configuration
In this study, four methods for simulating FD-OCT are compared: (1) original full-spectrum simulation with Mie functions [21] (referred as “Mie”); (2) simulation based on the conventional HG function and the corresponding IS method [25,38] (referred as “HG-IS”); (3) simulation with full-spectrum TT functions without IS (referred as “TT”); (4) our proposed method, denoted as “TT-IS”. Similar to the setting in [21], we also use water solution of polystyrene microspheres with the background refractive index (RI) ${n_m} = 1.33$ and particle RI ${n_s} = 1.58$ to form the simulated phantoms. The common configurations on phantom models, photon propagation, and optical parameters are summarized in Table 1.
Table 1. Common configuration details of simulations. N/A means “not available” or “not applicable”.
The scattering coefficients for each phantom are determined by Mie theory with particle size, density, RI, and wavelength as variables [13]. The probability of the additional scattering bias is set to 0.5 for methods with IS [25]. As the HG-IS method does not incorporate the wavelength dependency, the scattering coefficient is set corresponding to the commonly used wavelength of 1310 nm. Additionally, the anisotropy factor $g$ for HG-IS is set to 0.9 [39]. To quantitatively evaluate the performance of the following studies, we use the relative mean square error (rMSE) [40] $\hat e$ to measure the relative energy loss. Mathematically, the rMSE $\hat e$ for the 2-D data is defined as,
The whole program is run on a high-performance computing cluster equipped with Intel Xeon Gold 6258R Processor (2.7 GHz, 56 cores).
4.2 Fitting accuracy between the TT function and Mie function
We first investigate the ability of the five-parameter TT function to accurately capture the wavelength dependency in the Mie function. When employing NSLM to fit the TT function of unknown parameters and Mie function with Eq. (11), we constrain the values of each parameter in the TT function within the ranges specified in Eq. (1), and $A_{\rm max}=5$ in Eq. (13). To avoid constantly increasing the value of $\alpha _\cdot$ in the fitting algorithm, we set the range of $\alpha _f$ and $\alpha _b$ to $-{1}/{2}<\alpha _\cdot \leq \alpha _{\rm max}, \alpha _{\rm max}=10$. The initial values of ${\textbf {U}}_\lambda$ at each $\lambda$ are set to {0.5, 0.5, 1, 1, 0.99}.
The accuracy of the fitting is evaluated by Eq. (16), where $f(\cdot )$ is the TT function map with fitted parameters swept along $\lambda$ and $\theta$, and $f_{\rm GT}(\cdot )$ is the corresponding Mie function map. The number of sampling points of the scattering angle and wavelength is set to $N=3601$ and $M=1024$, respectively.
4.3 Validation: phantom simulation
To further validate our IS-accelerated method, we conduct some phantom-based simulations. The overview of A-line and B-scan simulations of four methods is shown in Table 2. The fitted reflectance map $R(l,\lambda )$ defined in [21] will be utilized as the representation of F-BDR, which is swept over the wavelength $\lambda$ and optical path length (OPL) $l$. For methods without IS, the F-BDR map is computed as,
where $\Delta l$ is the grid interval [21], and $W_\lambda ({l^*})$ is the whole energy of collected photons at OPL $l^*$ for wavelength $\lambda$. For TT-IS method, the F-BDR map will be modified as,Table 2. Summary of phantom simulations of FD-OCT. “A-mode” means the A-line simulation in FD-OCT, while “B-mode” is the simulation of B-scan; N/A indicates that this method is not suitable for this simulation or the computational cost of the simulation is too high.
4.3.1 A-mode phantom simulation
We conduct A-mode simulations on two types of phantoms: two semi-infinite phantoms (Phantom1 and Phantom2) and the multi-layer phantom, which is a stack of two phantoms of different thicknesses. As the HG-IS method does not incorporate the wavelength dependency, this method was not implemented.
For Mie and TT methods, we launch 20 million photon packets per wavelength for both simulation on the semi-infinite phantoms and the multi-layer one, which adds up to 20.48 billion for each simulation of F-BDR. Since the detection efficiency will be significantly increased by the importance sampling, we only use 5 thousand photon packets per wavelength for our proposed method (TT-IS) to simulate the semi-infinite phantoms and 10 thousand per wavelength for the multi-layer phantom. The quantitative results are measured by Eq. (16), in which $f(\cdot )$ is the F-BDR map, N = 1024 and M = 1024.
4.3.2 B-mode phantom simulation
We simulate 512 A-lines spanning a transverse range of [-1 mm, 1 mm] to generate a B-scan. Due to the excessive computational cost of collecting enough backscattered photons in methods without IS (Mie and TT), we had to abandon simulations using these two methods. Therefore, we focused on simulating the F-BDR and phantom images using the TT-IS method, while providing the image simulation with the HG-IS method as a comparison. For each A-line, 5 thousand photon packets per wavelength will be launched for both methods.
4.3.3 Volumetric phantom simulation
To further demonstrate the scalability of our proposed method, a volume of 128 frames (B-scans) is simulated with TT-IS method. For each B-scan, we simulate 512 A-lines within [-1 mm, 1 mm] and for each A-line, 2.5 thousand photon packets per wavelength will be launched.
5. Results
5.1 Fitting accuracy between the TT function and Mie function
Figure 3 presents the qualitative results of the fitting for two phantoms. To obtain a better view, the Mie and TT functions for each wavelength were normalized by the minimum value of all angles and then subjected to a logarithmic operation (see Fig. 3(a) and (b)). We also extract profiles at different wavelengths for both SPFs for a closer examination of the fitting accuracy, which are shown in polar coordinates in Fig. 3(c) and (d). The results demonstrate that by adjusting the parameters of the TT function, the trend of the Mie SPF at each wavelength can be well-fitted. However, the frequent fluctuations of the Mie function cannot be well fitted, since the TT function is a smooth analytic function in a closed form. For the fitting on Phantom1, as shown in Fig. 3(a), the fitted rMSE was 2.4%, whereas for Phantom2 (see Fig. 3(b)), the rMSE was 5.4%, slightly larger than the value for Phantom1. This may be attributed to more fluctuations in the Mie SPF due to the increase in particle radius.
Fig. 3. Results of the fitting of TT functions and Mie functions. (a) and (b) shows the Mie function map and the fitted results of TT from Phantom1 and Phantom2, respectively. For better comparison, profiles at four wavelengths are extracted in (c) and (d) for two phantoms. It is evident that TT function fits the trend of the Mie function well in wavelength and angle scanning. I - IV: 1270 nm, 1290 nm, 1310 nm, 1330 nm.
5.2 Validation: phantom simulation
5.2.1 A-mode phantom simulation
The normalized F-BDR maps of the semi-infinite Phantom1 and Phantom2 are shown in Fig. 4 and Fig. 5, respectively. Simulation results of the multi-layer phantom are presented on a logarithmic scale in Fig. 6. The extracted wavelength-resolved reflectance (top panels), as well as the depth-resolved profiles (right panels), facilitates a clear comparison of the results obtained using the Mie, TT, and TT-IS methods (corresponding to (b), (c), (d) of each figure).
Fig. 4. Simulation results of an A-line of Phantom1. (a) The implementation of the infinite Phantom1; (b), (c) and (d) are the normalized F-BDR maps of Mie, TT and TT-IS methods, respectively. The profiles at OPL of 250 $\mu {\rm m}$ (black) and 500 $\mu {\rm m}$ (red) are extracted to show the wavelength-dependent variation (top panels), while the OPL-resolved reflectance at wavelengths of 1270 nm (blue), 1310 nm (black) and 1330 nm (red) are also extracted for a closer look on the attenuation (right panels).
Fig. 5. Simulation results of an A-line of Phantom2. (a) The implementation of the infinite Phantom2; (b), (c) and (d) are the normalized F-BDR maps of Mie, TT and TT-IS methods, respectively. Wavelength-resolved and OPL-resolved profiles are extracted in the same setting of Fig. 4.
Fig. 6. Simulation results of an A-line of the multi-layer phantom. (a) The implementation of the multi-layer phantom; (b), (c) and (d) are the normalized F-BDR maps of Mie, TT and TT-IS methods shown on the logarithmic scale, respectively. Different from the Fig. 4 and Fig. 5, the profiles at OPL of 500 $\mu {\rm m}$ (red) and 1000 $\mu {\rm m}$ (black) are extracted to highlight the distinct variation of two types of layers along the wavelength (top panels). The results obtained by TT-IS (d) are consistent with the F-BDR variation and the structural features of different layers obtained through the Mie method (a).
Qualitatively, the different wavelength dependency exhibited by the two phantoms is well preserved in both the TT method and the proposed TT-IS method, which is very close to the ground truth from the wavelength direction. Moreover, a comparison of the TT-IS method (Fig. 6(d)) with the Mie method (Fig. 6(b)) reveals that the structure of the multiple layers and the positions of boundaries are correctly visible. However, the TT and TT-IS methods seem to slightly overestimate the BDR intensity in the deep region.
Table 3 presents the calculation time and rMSE of the three simulations. Quantitatively, although the TT-IS method tends to overestimate the strength of F-BDR in the deep region, its energy error relative to the Mie method is less than 5% in all three simulations, especially less than 2% on the multi-layer phantom that launches more photons. Remarkably, the TT-IS method achieves a speedup of 573$\times$ for Phantom1, 194$\times$ for Phantom2, and 373$\times$ for the multilayer phantom, compared to the Mie method.
Table 3. Quantitative results of the A-line simulations. The TT-IS method achieves a remarkable acceleration of over two orders of magnitude compared to the ground truth Mie method, with the required time highlighted in bold.
5.2.2 B-mode phantom simulation
Thanks to the improved simulation speed, the calculation of a B-scan consisting of 512 A-lines can now be completed in a vastly improved time. The F-BDR simulation of this B-scan is shown in Fig. 7. We simulate 512 A-lines in parallel, each taking an average of only 0.4 hours. The total simulation time adds up to 204.8 CPU core hours. Figure 7(b) displays the F-BDR of a phantom with the structure depicted in Fig. 7(a), which is formed as a three-dimensional data cube along the wavelength, OPL, and transverse positions. It demonstrates the ability to resolve different layers and spherical structures, with clear variations of BDR observed along the OPL direction. To better demonstrate the preserved wavelength dependency, a threshold was applied to the data cube to only retain the upper surface portion of the spherical structure, as indicated by the red dashed box in Fig. 7(b), as presented in Fig. 7(c). The varying distribution of intensity depicted in Fig. 7(c) showcases the wavelength dependency and the trend observed is similar to profiles presented in the top panel of Fig. 5(d).
Fig. 7. F-BDR simulation results of a 2-D structured phantom. (a) The structure of the simulated phantom. 512 A-lines are calculated within a transverse range of [-1 mm, 1 mm]; (b) Simulated F-BDR of the B-scan, which is formed as a 3-D data cube swept along the wavelength, OPL, and transverse position; (c) F-BDR of the upper surface of the spherical structure after the thresholding. The wavelength dependency of the F-BDR can be clearly visible.
Using the obtained F-BDR, we further reconstruct the B-mode image $I(x,z)$ of FD-OCT, as shown in Fig. 8(b), where $x$ represents the transverse position and $z$ represents the axial depth, which is half of the OPL. Additionally, the image obtained with the HG-IS method, which is similar to [38], is presented in Fig. 8(a). It is evident that, compared to the method based on the HG function and its corresponding IS, the proposed method better displays the structural characteristics of the phantom and produces more realistic speckle patterns.
Fig. 8. Reconstructed B-scan image of the phantom implemented in Fig. 7(a). (a) Image obtained from the HG-IS method; (b) Image reconstructed with the F-BDR cube obtained from the proposed TT-IS method, which is displayed in Fig. 7(b).
5.2.3 Volumetric phantom simulation
A complex 3-D phantom is built as shown in Fig. 9(a). The layered structure is implemented using Phantom1, while a branched structure is constructed with Phantom2. The empty space is the background medium (water). Figure 9(b) shows the simulated volume of the phantom. The overall simulation time is 5898 CPU core hours. Our method has simulated a realistic volume, and some phenomena, such as shadows below the strongly scattered medium (Phantom2) can be observed. It should be noted that the full-spectrum volume simulation is almost impossible with the original simulator without IS-based acceleration.
Fig. 9. (a) 3D Phantom implementation containing 128 frames; (b) Simulated FD-OCT volume of the phantom.
6. Discussion
6.1 Error propagation of the algorithm
The proposed IS-accelerated simulation method is divided into two steps. The first step is to fit Mie functions and TT functions with unknown parameters, and the second step is to apply the derived TT-B function to perform biased sampling in the scattering process. Each step brings errors.
The first step of the proposed method can introduce errors mainly due to two factors: (1) the inaccuracy of the fitting procedure, and (2) the loss of a large number of fluctuations in the Mie function. Both could lead to a slight overestimation in the fitting values of the backscattering part in TT SPF, which results in the misestimation of BDR intensity in the deep region. The former can be mitigated by adjusting the weighting method specified by $h(\cdot )$ in Eq. (11) or by using a more robust fitting method than NLSM, such as genetic algorithms [41]. For the latter, due to the smoothness of the TT function, the fitting loss of fluctuation in the Mie function seems inevitable, and the specific effects of the lack of these fluctuations on scattering need to be further studied.
In the second step, the errors introduced by the TT method are further amplified by the TT-IS method. For example, by comparing (b) and (c) as well as (b) and (d) in Fig. 4 and Fig. 5, it can be observed that when the TT method results in certain F-BDR peak shift and energy overestimation, TT-IS will further exacerbate these phenomena, thereby introducing more errors. This is due to the improved sampling efficiency of TT-IS, as well as the fact that the five parameters of TT are simply copied into TT-IS. Adjusting the parameters in TT-IS may help to reduce these errors, but in cases where little energy loss occurs, such adjustments may be unnecessary.
6.2 Difference in the B-scan image between HG-IS and TT-IS method
The B-scan images presented in Fig. 8 exhibit a significant difference. Our proposed method produces clear structures of the phantom, while the HG-based method fails to do so. It is because of the misestimation of the backscattering component in the HG function, as well as the fact that the HG function does not account for wavelength dependency and employs the same anisotropy factors across different phantoms. Although the scattering coefficients of the two types of phantoms utilized in the HG-based method differ, it appears that SPF plays a more significant role in the scattering process, resulting in an unclear central spherical structure in the HG-based image (see Fig. 8(a)).
6.3 Generalization in BDR-based applications
Although our proposed method has only been validated for use in FD-OCT, it can be easily adapted to other BDR-dependent applications. The proposed TT-B function for controlling biased scattering is based on the TT function, which exhibits excellent generalization properties and can be applied in other areas without imposing additional conditions. The proposed biased scattering process is adapted from a previous study [25] that directs photons towards the detector, resulting in a biased distribution. However, we believe that for different applications, the biasing process can be modified to direct photons towards specific regions of interest. In such cases, the TT and TT-B functions can still be used for unbiasing and biasing scattering sampling, respectively.
6.4 Further speedup with CUDA parallel computing
Although the simulation combined with IS acceleration is faster than the original simulator [21], it still takes hundreds to thousands of core hours to produce a well-simulated image or volume. Considering the potential applications of the MC simulation, such as inverse estimation of optical parameters [42,43], requires lower computational cost, we need to further improve the simulation speed. Using Compute Unified Device Architecture (CUDA) for parallel computing can accelerate the MC simulation hundreds of times while ensuring accuracy [26]. In the future, we will migrate the proposed algorithm to the CUDA platform for further acceleration.
6.5 Application on biological tissues
While the proposed method exhibits promising performance in simulating particle-based phantoms, validating it with real biological tissues poses potential challenges. One challenge is that biological tissues are based on complex refractive index fluctuations rather than distributions of spherical particles. In such cases, the smooth TT SPF may offer a more effective simulation of the scattering properties of biological tissues compared to the Mie SPF. Exploring methods to establish connections between particle-based models and complex refractive index fluctuations would be one of our future research topics. Secondly, it is difficult to do the parameter fitting of the TT SPF with the Mie SPF due to the unknown particle distribution. In this case, it may be necessary to experimentally measure the TT SPF.
6.6 Progressive acceleration of the algorithm
From Table 3, it is apparent that both TT and TT-IS methods achieve the acceleration of the computational time compared to the Mie method. Since the TT SPF is smoother than Mie SPF, it is easier to inversely sample the scattering angle using the TT function, which makes the TT method a little faster than the Mie method. However, this acceleration is also limited by the collection efficiency of backscattered photons. Only by adding IS to the TT method can we improve the collection efficiency and greatly reduce the computational cost.
7. Conclusion
To tackle the problem of the excessive computational cost of MC simulation for F-BDR, we propose a novel IS-based acceleration scheme utilizing the TT function, which can accurately approximate the Mie SPF. We derive a TT-based general biased SPF, TT-B function, tailored for the case of biasing photons towards the region of interest (detector). By fitting the TT function and Mie function, we obtain wavelength-dependent parameters for both the TT and TT-B functions. These functions are then integrated into the IS process to bias photons towards the detector, resulting in an efficient and accurate simulation of backscattered events at each wavelength. The application on the full-spectrum simulator of FD-OCT fully demonstrates the effectiveness of the proposed IS-accelerated F-BDR simulation method: compared with the original simulator, the F-BDR of a multi-layer phantom can be produced in 373$\times$ speed with the rMSE less than only 2%. Moreover, this method can also be extended to the simulation to cross-sections and volumes, which makes it a powerful tool for studying numerous forward modeling and inverse problems in various applications that require F-BDR simulation.
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. The reviewers’ valuable comments were much appreciated, which greatly improved the manuscript.
Disclosures
The authors declare no conflicts of interest.
Data availability
The source codes of this manuscript are publicly available online at [33].
References
1. W. Drexler and J. G. Fujimoto, Optical Coherence Tomography: Technology and Applications (Springer Science & Business Media, 2008).
2. C. S. Mulvey, C. A. Sherwood, and I. J. Bigio, “Wavelength-dependent backscattering measurements for quantitative real-time monitoring of apoptosis in living cells,” J. Biomed. Opt. 14(6), 064013 (2009). [CrossRef]
3. N. Bosschaart, D. J. Faber, T. G. van Leeuwen, and M. C. Aalders, “Measurements of wavelength dependent scattering and backscattering coefficients by low-coherence spectroscopy,” J. Biomed. Opt. 16(3), 030503 (2011). [CrossRef]
4. 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]
5. J. Shin, B. T. Bosworth, and M. A. Foster, “Single-pixel imaging using compressed sensing and wavelength-dependent scattering,” Opt. Lett. 41(5), 886–889 (2016). [CrossRef]
6. V. V. Tuchin, “Tissue optics and photonics: light-tissue interaction,” J. Biomed. Photonics Eng. 1, 98–134 (2015). [CrossRef]
7. C. Photiou, E. Bousi, I. Zouvani, and C. Pitris, “Using speckle to measure tissue dispersion in optical coherence tomography,” Biomed. Opt. Express 8(5), 2528–2535 (2017). [CrossRef]
8. M. Rinehart, Y. Zhu, and A. Wax, “Quantitative phase spectroscopy,” Biomed. Opt. Express 3(5), 958–965 (2012). [CrossRef]
9. 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]
10. 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]
11. C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt. 18(5), 050902 (2013). [CrossRef]
12. 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]
13. H. C. Hulst and H. C. van de Hulst, Light Scattering by Small Particles (Courier Corporation, 1981).
14. X. Wang and L. V. Wang, “Propagation of polarized light in birefringent turbid media: a Monte Carlo study,” J. Biomed. Opt. 7(3), 279–290 (2002). [CrossRef]
15. C. He, H. He, J. Chang, B. Chen, H. Ma, and M. J. Booth, “Polarisation optics for biomedical and clinical applications: a review,” Light: Sci. Appl. 10(1), 194 (2021). [CrossRef]
16. R. Hennessy, S. L. Lim, M. K. Markey, and J. W. Tunnell, “Monte Carlo lookup table-based inverse model for extracting optical properties from tissue-simulating phantoms using diffuse reflectance spectroscopy,” J. Biomed. Opt. 18(3), 037003 (2013). [CrossRef]
17. P. Naglič, F. Pernuš, B. Likar, and M. Bürmen, “Limitations of the commonly used simplified laterally uniform optical fiber probe-tissue interface in Monte Carlo simulations of diffuse reflectance,” Biomed. Opt. Express 6(10), 3973–3988 (2015). [CrossRef]
18. G. M. Palmer and N. Ramanujam, “Monte-Carlo-based model for the extraction of intrinsic fluorescence from turbid media,” J. Biomed. Opt. 13(2), 024017 (2008). [CrossRef]
19. C. Liu, N. Rajaram, K. Vishwanath, T. Jiang, G. M. Palmer, and N. Ramanujam, “Experimental validation of an inverse fluorescence Monte Carlo model to extract concentrations of metabolically relevant fluorophores from turbid phantoms and a murine tumor model,” J. Biomed. Opt. 17(7), 0780031 (2012). [CrossRef]
20. 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]
21. J. Mao, Y. Ling, P. Xue, and Y. Su, “Monte Carlo-based full-wavelength simulator of Fourier-domain optical coherence tomography,” Biomed. Opt. Express 13(12), 6317–6334 (2022). [CrossRef]
22. J. Schmitt and K. Ben-Letaief, “Efficient Monte Carlo simulation of confocal microscopy in biological tissue,” J. Opt. Soc. Am. A 13(5), 952–961 (1996). [CrossRef]
23. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef]
24. 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]
25. 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]
26. 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(08), 083012 (2022). [CrossRef]
27. R. Y. Rubinstein and D. P. Kroese, Simulation and the Monte Carlo Method (John Wiley & Sons, 2016).
28. R. Yuan, J. Ma, P. Su, Y. Dong, and J. Cheng, “Monte-Carlo integration models for multiple scattering based optical wireless communication,” IEEE Trans. Commun. 68(1), 334–348 (2020). [CrossRef]
29. L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” The Astrophys. J. 93, 70–83 (1941). [CrossRef]
30. S. L. Jacques, C. Alter, and S. A. Prahl, “Angular dependence of HeNe laser light scattering by human dermis,” Lasers Life Sci 1(4), 309–333 (1987).
31. M. Canpolat and J. R. Mourant, “High-angle scattering events strongly affect light collection in clinically relevant measurement geometries for light transport through tissue,” Phys. Med. Biol. 45(5), 1127–1140 (2000). [CrossRef]
32. S. L. Jacques and N. J. McCormick, “Two-term scattering phase function for photon transport to model subdiffuse reflectance in superficial tissues,” Biomed. Opt. Express 14(2), 751–770 (2023). [CrossRef]
33. J. Mao, “fullwaveoct,” Github, 2023, https://github.com/Jianing-Mao/fullwaveOCT.
34. L. O. Reynolds and N. J. McCormick, “Approximate two-parameter phase function for light scattering,” J. Opt. Soc. Am. 70(10), 1206–1212 (1980). [CrossRef]
35. L. Wang, S. L. Jacques, and L. Zheng, “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef]
36. L. V. Wang and H.-i. Wu, Biomedical Optics: Principles and Imaging (John Wiley & Sons, 2012).
37. P. E. Gill and W. Murray, “Algorithms for the solution of the nonlinear least-squares problem,” SIAM J. Numer. Anal. 15(5), 977–992 (1978). [CrossRef]
38. S. Zhao, “Advanced Monte Carlo simulation and machine learning for frequency domain optical coherence tomography,” Ph.D. thesis, California Institute of Technology (2016).
39. S. L. Jacques and L. Wang, “Monte Carlo modeling of light transport in tissues,” Optical-Thermal Response of Laser-Irradiated Tissue pp. 73–100 (1995).
40. E. Holst and P. Thyregod, “A statistical test for the mean squared error,” J. Stat. Comput. Simul. 63(4), 321–347 (1999). [CrossRef]
41. J. H. Holland, “Genetic algorithms,” Sci. Am. 267(1), 66–72 (1992). [CrossRef]
42. G. M. Palmer and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. part i: Theory and validation on synthetic phantoms,” Appl. Opt. 45(5), 1062–1071 (2006). [CrossRef]
43. K. C. Chong and M. Pramanik, “Physics-guided neural network for tissue optical properties estimation,” Biomed. Opt. Express 14(6), 2576–2590 (2023). [CrossRef]
