1. Introduction

Optical coherence tomography (OCT) is a high-resolution, non-invasive imaging technique increasingly used for medical and biological applications. It is based on the principle of low-coherence interferometry, which relies on utilizing a light source with a broad wavelength bandwidth to ensure coherence gating, thus leading to a high axial resolution of a few microns [1]. In OCT systems, an interference spectrum is created by combining the reflected light from a reference mirror with that from the sample arm. This pattern is used to obtain a reflectance profile of the sample as a function of depth, namely an OCT A-line. OCT imaging primarily relies on detecting ballistic and quasi-ballistic backscattered photons from turbid environments, such as biological tissues. As light penetrates deeper into the tissue, its intensity decreases exponentially due to scattering and absorption, subsequently limiting the imaging depth to a few millimeters. OCT is now a standard diagnostic tool in the field of ophthalmology and a promising candidate in various other fields, including cardiology, dermatology, gastroenterology, and laryngology [2–5].

Light propagation in a turbid medium can be modeled in simple geometries by analytically solving the radiative transfer equation (RTE) [6]. However, when imaging biological tissue with high scattering properties, the mean free path of light within the tissue is short, which renders the RTE calculation and its inverse problem reconstruction too complex. In such cases, numerical simulation techniques using Monte Carlo (MC) methods have emerged as the gold standard for modeling photon interactions in biological tissue [7–12]. With MC methods, light is considered a collection of photon packets tracked as they undergo multiple events of reflection, refraction, absorption, and scattering while traveling through the turbid medium. Over the past years, several MC algorithms have been developed, amongst which Monte Carlo multi-layered (MCML) is the most widely known, forming the foundation of numerous other simulators for biomedical applications [9]. In MCML, an ensemble of photon packets is launched as a pencil beam and perpendicularly to the surface of a multi-layered medium, with a simulated beam geometry that is cylindrically symmetric around the $z$-axis. The use of voxels, tetrahedrons, and 3D volumetric meshes further allowed the extension of planar and cylindrical geometries to model more realistic tissue structures [11,13–16]. To accelerate the speed and improve computational efficiency, it was also proposed to use a graphics processing unit (GPU), enabling simultaneous simulations of a higher number of photon packets [17,18].

Monte Carlo methods were previously applied to simulate OCT signals [10,19–26] and some open-source OCT simulators recently became available [27–29]. Monte Carlo simulation for tissue layers with embedded objects for OCT (MCEO-OCT) is a modified version of MCML that embeds geometric shapes in the multi-layered tissue to mimic vessels and utilizes the method of importance sampling to accelerate the computation time [27]. The OCT massively parallel simulator (OCT-MPS) is based on an MC approach that exploits a tetrahedron mesh and GPU programming [28]. These two simulators were developed for time-domain OCT applications as only photons describing a path length difference with the reference path length smaller than the coherence length of the OCT system are tagged as contributing to the signal. For Fourier-domain applications, three MC simulators have previously been developed: one uses a layered tissue geometry with Gaussian beam focusing [30], another exploits a 3D voxel system with a pencil beam [31], while the third tool proposed an acceleration scheme based on importance sampling with a phase-function that depends on the wavelength [29]. The latter tool is open-source. These tools can be used to understand OCT signals better, aid in the instrumental design and development of OCT-based applications, guide the development of OCT signal processing algorithms, and more recently train neural networks to improve OCT image reconstruction [10,19–25,27,28,30–33].

In this work, we propose an open-source Fourier-domain simulator, MCOCT, modified from MCML 3D voxel geometry mcxyz, in which we implemented modifications to accommodate the particular requirements of Fourier-domain OCT. To achieve depth-resolved reflectometry in addition to simulating photon propagation, photon packets are simulated at a single central wavelength and tracked to record their traveled distance and extract phase as a function of wavelength. This strategy allows mimicking coherence gating. A photon detection system with an angle of acceptance was also implemented to mimic the numerical aperture of the imaging system [10]. To achieve a significant signal-to-noise ratio and expedite computation time, importance sampling was also employed [24,34]. Gaussian illumination was further incorporated to model the light intensity profile and better represent OCT imaging. For numerical validation, the optical fluence rate generated by MCOCT was compared to the one simulated by Monte Carlo eXtreme (MCX) [18] using the same optical properties and geometry. For experimental validation, OCT signals generated by MCOCT were compared to OCT signals acquired at specific numerical apertures in synthetic phantoms with known optical properties.

2. Material and methods

In the following subsections, we will sequentially describe the specific parts of the MCOCT simulator.

2.1 Monte Carlo photon propagation

A standard MC simulation starts with the launch of a photon packet into a medium that is defined by its index of refraction $n$, absorption coefficient $\mu _ {a}$, scattering coefficient $\mu _ {s}$, and anisotropy coefficient $g$ [9]. For voxel-based approaches, each voxel is associated with a tissue type and its optical properties. Photon packets randomly propagating through the geometry are associated upon launch with a weight $W=1$ representing the packet intensity that decreases as it undergoes events. At the site of an absorption interaction, $W$ is updated according to $W = W \cdot \exp {-(\mu _a\cdot s)}$, where $s$ is the distance traveled by the photon in that step [9]. In the case of a scattering interaction, the diffusion direction is calculated from the azimuthal angle $\phi$ and the deflection angle $\theta _ {s}$ [6]. $\phi$ is uniformly distributed in the range $[0,2\pi ]$, whereas $\theta _ {s}$ is typically modeled by the Henyey-Greenstein (HG) function $f_{HG}$ that depends on the optical properties of the medium [35,36]. Since most biological tissues have an anisotropy coefficient $g \in [0.90, 0.95]$ at the clinically applied wavelengths of OCT, the photons are more prone to scatter in the forward direction [37]. When $W$ is below a predetermined critical threshold, a roulette process begins in which the photon packet is associated with a probability to continue its propagation, otherwise, the sequence is terminated for this particular photon packet, and $W$ is zeroed [9]. The simulation ends upon the termination of all photon packets.

2.2 OCT signal simulation

In order to simulate OCT A-lines, we introduced a new beam-focusing method based on Gaussian beams, used statistical methods to increase the probability of photon detection, tracked photons path lengths, and incorporated a system mimicking photon illumination and detection.

2.2.1 Gaussian beam description

In OCT, beam focussing impacts both signal intensity and resolution. The diffraction-limited focal spot of a Gaussian beam at a certain wavelength $\lambda$ is determined by its beam waist $w_0$. Employing the paraxial approximation, the beam waist can be described as $w_0 = 2\lambda f/(\pi \cdot D)$, relative to the focal length of the imaging lens $f$, the beam center wavelength $\lambda$ and the beam diameter at the lens $D$ [38]. In OCT imaging, light propagates through several layers of different indices of refraction $n_i$ with each layer affecting the axial position and the width of the beam waist, i.e., the beam curvature. Figure 1 illustrates the curved trajectory of such a beam depicting the individual parameters where $z=0$ corresponds to the sample surface, described by the red line. The position of the focus in the tissue is given by $z_f$ and the distance between the tissue surface and imaging lens by $z_{f,I}$, presenting a more general approach compared to assuming the focusing lens to be at $z=0$ [38].

 

Fig. 1. Description of the geometry and parameters used to simulate the Gaussian focusing beam. The dashed lines show the trajectory in a crossing propagation. Variables are described in Table 5 of the Appendix.

Download Full Size | PDF

To implement this beam curvature in the MC method, the beam was focused by propagating photons along the trajectory by a small step unit $s_\text {min}$ to the next position $(x,y,z)$. Thereby, the unit of the step size is defined by the optical properties of the tissue. After each step, the direction unit vector $\hat {u}$ is updated to ensure precise photon tracking along the initial trajectory. Considering a photon propagation through two layers of differing refractive indices $n_1$ and $n_2$, the beam waist of the initial beam in a medium with $n_2$ can be described as

The Rayleigh range of the beam $z_R$ in the medium of $n_2$ is

The photon packet propagates by $s_{\textrm{min}}$ until it undergoes its first scattering event, after which it propagates in accordance with a traditional MC simulation [9]. To allow the photon packet to follow a curved path, the packet is moved by a distance $s_{\textrm{min}}$ before its direction $\hat {u}$ is recalculated using Eq. (4). The photon is moved this way until it has traveled the entire distance $s$. To assess the performance of the Gaussian illumination, simulations were performed with a pencil beam and a geometric focused beam [9,18].

2.2.2 Importance sampling and photon splitting methods

To increase the probability of backscattering events, the method of importance sampling was implemented and combined with photon splitting (see flowchart in Fig. 2) [24,34]. Importance sampling is a probabilistic method to accelerate MC simulations for OCT by biasing random events leading to backscattering, thereby increasing the probability of photon detection and thus enhancing computational efficiency. An updated biased scattering function $f_B$ based on the HG probability density function $f_{HG}$ was then introduced. After a scattering event, this function modified the photon propagation direction to be oriented towards the detector, i.e., the collecting fiber system. The resulting bias was corrected by introducing a corresponding likelihood ratio as

 

Fig. 2. Flowchart describing the scattering process of the simulator. The randomly generated number is noted $\xi _i\in [0,1]$ while $p$ is the probability for biased versus unbiased scattering.

Download Full Size | PDF

2.2.3 Optical configuration and tracking algorithm

Optical illumination and detection configurations were simulated using illumination and photon collection modules with identical properties of numerical aperture (NA), geometry, and spatial locations. Simulations of multiple A-lines were achieved by varying the lateral positions of the illumination and detection fibers. When the backscattered photon returned within the fiber’s NA and radius, the photon was labeled as detected. The traveled distance of the photon packet in each different layer was tracked in addition to the fluence absorbed by the sample. This strategy allowed us to adjust the phase of the signal for media with different refractive indices afterward. Upon photon detection, the cumulative traveled distance (total pathlength $TP$), weight $W$, and likelihood $L$ of each photon were saved.

This method was implemented by modifying an existing open-source C-language software package that is available on the website of the Oregon Medical Laser Center [13]. The MCOCT photon tracking algorithm flowchart is shown in Fig. 3 and describes the steps of a standard MC simulation (gray background) as well as our methodological contributions (white background).

 

Fig. 3. MCOCT flowchart. Variables are described in Table 5 of the Appendix.

Download Full Size | PDF

2.3 OCT signal reconstruction

For each A-line defined by a lateral ($x$) position of the illumination/detection fiber module, a number $N$ of detected photons with their corresponding weight, likelihood, and total pathlength was acquired. These parameters were then used to reconstruct the OCT signal following a series of mathematical processes described in Fig. 4. To calculate the signal across the spectral bandwidth $\Delta \lambda$ at a central wavelength $\lambda _0$, a discrete $k$ array corresponding to the $m^{\text {th}}$ wavenumber ($k_m = 2\pi n/\lambda _m$) was created with $\lambda _\text {min} = \lambda _0 - \frac {\Delta \lambda }{2}$ and $\lambda _\text {max} = \lambda _0 + \frac {\Delta \lambda }{2}$. The collected signal $S_S(k_m)$ from the sample was then reconstructed as a function of wavelength and for all $j$ detected photon packets using this expression:

The reference arm signal amplitude $S_R(k_m)$ was calibrated as the average of the modulus $|\cdot |$ of the signal amplitude of the sample arm to ensure comparable power in both arms :

By emulating a balanced detection configuration, the signal at the detector was generated by interfering the sample arm signals $S_S(k_m)$ with the reference arm signals $S_R(k_m)$ and therefore creating the interference fringe signal $I(k_m)$ defined by

Consequently, $I(k_m)$ was defined as a function of the wavelength and required further Fourier-domain processing to generate the OCT signal as a function of depth by applying a Hanning window, zero padding, and discrete Fourier transform. These processes were performed in Matlab (Version R2021a, Natick, MA, USA).

 

Fig. 4. Optical coherence tomography (OCT) signal reconstruction. Definitions of variables: $TP_N$, total pathlength; $W_N$, weight; $L_N$, likelihood; $N$, $N^{\text {th}}$ collected photon packets; $k$, wavelength; $S_S(k_m)$, signal amplitude from the sample arm; $S_R(k_m)$, signal amplitude from the reference arm; $I(k_m)$, reconstructed OCT signal; $k_m$, $m^{\text{th}}$ wavenumber.

Download Full Size | PDF

2.4 Tissue-mimicking optical phantoms

Optical phantoms based on the optical properties of biological tissues were fabricated and characterized to validate this simulator.

2.4.1 Fabrication

To validate our simulator, optical phantoms mimicking biological tissue were developed with polydimethylsiloxane (PDMS) matrices. PDMS was selected for an accurate representation of the refractive index $n\sim 1.4$ of biological tissue [39]. The matrix demonstrates good durability and is visible with the OCT system at a wavelength of $1310$ nm. Nigrosine was utilized as the absorbing agent and titanium dioxide powder as the scattering agent to mimic biological tissue (see Table 1). Three one-layer phantoms and one bi-layer phantom were fabricated with different optical properties by varying quantities of nigrosine and titanium dioxide. The bottom layer of the bi-layer phantom was created using the same recipe as the one-layer phantom #3, wherefore the properties of these layers are assumed equal.

Table 1. Quantities of the scattering (titanium dioxide) and absorbing (nigrosine) agents of the phantoms developed with polydimethylsiloxane (PDMS) matrices.

View Table | View all tables in this article

2.4.2 Optical properties characterization

The absorption coefficient $\mu _a$ and reduced scattering coefficient $\mu _s'$ defined as $\mu _s' = (1-g) \mu _s$ were determined for the one-layer phantoms. To acquire the transmission and reflection spectra, a LAMBDA 1050 spectrophotometer from PerkinElmer (Waltham, MA, USA) was used. This commercial device is equipped with a single integration sphere capable of receiving both diffuse and specular lights. The collected spectra were processed using the inverse adding-doubling method to estimate $\mu _a$ and $\mu _s'$ [40,41]. The determination of the anisotropy coefficient $g$ was unattainable through direct measurement and was thus estimated by fitting OCT A-lines measured experimentally (see Sec. 2.4.3) to a theoretical model. In a first approximation, assuming a homogeneous turbid medium, the signal amplitude $I(z)$ at a depth $z$ can be modeled as a single exponential decay function following a variant of the Beer-Lambert law [42–44]:

Table 2. Optical properties of the phantoms characterized by spectrophotometry: $\mu _a$, absorption coefficient; $\mu _s'$, reduced scattering coefficient; $n$, refractive index; $g$, anisotropy coefficient.

View Table | View all tables in this article

2.4.3 OCT imaging of phantoms

OCT imaging was performed with a VCSEL Fourier-domain system (Thorlabs, NJ, USA) at $\lambda _0 = 1310$ nm with $\Delta \lambda = 130$ nm. Two lenses of different focal lengths were used: a short one with $f = 30$ mm and a long one with $f=75$ mm. The position of the focus inside each optical phantom $z_f$ was determined separately for each focal lens and summarized in Table 3. For each experiment, OCT imaging generated 300 A-lines which were averaged for comparison with MCOCT.

Table 3. Position of focus $z_f$ [mm] below the phantom surface for each focal lens $f$.

View Table | View all tables in this article

2.5 MCOCT simulations for comparison with MCX simulations

To validate the beam focusing performance, the flux output (optical fluence rate) generated by MCOCT was compared to the one from MCX [18]. The accumulated energy was stored within the 3D voxel matrix during the simulation. The 3D matrix with the deposited weight was normalized relative to the input power, resulting in the local absorption rate per watt (W) of the incident light in units of W/cm³/W_incident. Dividing this rate by the absorption coefficients of the corresponding voxels resulted in the local fluence rate $\phi$ in W/cm²/W_incident. The tissue was represented in both simulators by a $200 \times 200 \times 300$ voxel geometry with isotropic voxels of $0.001~$cm³. The distance between the imaging lens and tissue surface $z_{f,I}$ was set at 0.85 cm with a beam width at imaging lens $D=0.1$ cm. In MCOCT, the illumination beam had a Gaussian profile at the imaging lens. In MCX, the trajectory was simulated as a hyperboloid of one sheet and the photon packets propagated in straight lines along the double-ruled surface, thereby mimicking a Gaussian-distributed photon launch [18]. For both simulations, optical properties were $\mu _a=0.001$ cm^-1, $\mu _s=0.01$ cm^-1, $n=1$, and $g=0.5$. MCOCT and MCX simulations were executed using 10⁶ photon packets and were run on a desktop computer equipped with an Intel Xeon E5 CPU at 3.50GHz, 16GB RAM and an NVIDIA TITAN Xp card.

2.6 MCOCT simulations for comparison with OCT imaging in phantoms

Optical properties from Table 2 were used to generate OCT A-lines with MCOCT. Importance sampling was incorporated with a bias coefficient of $a=0.9$ and an additional bias probability of $p=0.8$, as defined in [24,34]. A minimum step $s_\textrm{{min}} = 10^{-4} \,\text {cm}$ was used for these simulations. A scheme of 1000 equidistant illumination and detector pairs along the $x$-axis from $x=-1$ cm to $x=1$ cm was used and resulted in simulating 1000 A-lines for a single B-scan. To minimize simulation time, a one-voxel slab of $2\times 2\times 0.5$ cm³ was used to model the one-layer phantoms. For comparison with the bi-layer phantom, the medium was modeled as a two-voxel slab of $2\times 2\times 1.8$ cm³. For all simulations, the incident light beam diameter at the focusing lens $D$ was set to 0.07 cm for imaging at $\lambda _0 = 1310$ nm. For OCT signal reconstruction, $\lambda _\text {min} = 1245$ nm and $\lambda _\text {max} = 1375$ nm. Optical properties were approximated to be constant over the OCT reconstruction spectrum. All simulations were executed using between 10⁸ and 10⁹ photon packets and were ran on the same desktop computer as above. Simulated A-lines are further averaged among detectors (1000) to reduce speckle contrast.

3. Results

Figure 5(a) and (b) show the optical fluence rates $\phi$ at $y = 0$ mm generated by MCOCT and MCX, respectively. In both simulations, the beam appears curved as expected for a Gaussian distribution. The simulation time was 60 seconds for MCOCT and one second using MCX. Figure 5 (c) shows a very good agreement between optical fluence rates simulated along the $x$-axis at the depth of the beam focus ($z=1.5$ mm) by MCOCT (red) and MCX (blue). Figure 5 (d) presents the optical fluence rates simulated along the $z$-axis at a fixed $x$-position ($x=0$ mm) by MCOCT (red) and MCX (blue). Both curves progress similarly around the beam waist but show more discrepancies with distance including a fastest decrease for MCX compared to MCOCT.

 

Fig. 5. Optical fluence rates $\phi$ [W/cm²/W_incident] comparison between (a) MCOCT and (b) MCX simulated with the same optical properties. (c) Comparison of optical fluence rates along $x$ at $z=1.5$ mm and $y=0$ mm generated by MCOCT (red) and MCX (blue). (d) Comparison of optical fluence rates along $z$ at $x=0$ mm and $y=0$ mm simulated with MCOCT (red) and MCX (blue).

Download Full Size | PDF

Figure 6 and Fig. 7 show the comparison between A-lines generated by MCOCT (blue) and those measured in the 4 experimental phantoms (red) using an imaging lens with a focal length of $f=30$ mm and $f=75$ mm, respectively. At $f=30$ mm, the attenuation of the MCOCT A-line follows a similar pattern to experimental data in the one-layer phantoms except in a region close to the surface of the sample where the signal intensity is lower in experimental data. In the bi-layer phantom (phantom #4, (d)), the change of optical properties from the superficial layer to the deep layer is marked by a transient change in signal intensity at a depth corresponding to the interface between layers (0.9 mm). A change in the attenuation slope is also observed with the change of optical properties between layers. For all comparisons, the accuracy of the simulated data also decreased with increasing depth as a result of lower signal intensity due to fewer collected photons. At $f=75$ mm, signal intensity patterns were similar to that of $f=30$ mm while signal intensity close to the surface presented a more accurate qualitative match. For both focal lengths, A-lines obtained by simulating a pencil beam and a geometric focused beam differed significantly from the A-line simulated with the Gaussian beam as well as from experimental OCT signals acquired in Phantom #1 (see Fig. 8 in Appendix).

 

Fig. 6. Experimentally obtained (red) vs. simulated (blue) OCT A-lines using an imaging lens of $f=30$ mm of (a) Phantom #1, (b) Phantom #2, (c) Phantom #3, and (d) Phantom #4.

Download Full Size | PDF

 

Fig. 7. Experimentally obtained (red) vs. simulated (blue) OCT A-lines using an imaging lens of $f=75$ mm of (a) Phantom #1, (b) Phantom #2, (c) Phantom #3, and (d) Phantom #4.

Download Full Size | PDF

Table 4 reports computation times for generating a B-scan of the phantoms presented in Fig. 6 and Fig. 7 along with the quantity of simulated and detected photons. Detected photons varied between $1.5\times 10^5$ to $2.9\times 10^6$ for a number of photons launched varying between $10^8$ to $10^9$, yielding a ratio of detected/launched between 0.0525 and 0.33%.

Table 4. Computation time and quantity of launched and detected photons for MCOCT simulations presented in Fig. 6 and Fig. 7.

View Table | View all tables in this article

4. Discussion

In this work, an OCT simulator (MCOCT) based on Monte Carlo methods was implemented and validated. The simulator was based on a Gaussian illumination profile and the importance sampling method as well as signal reconstruction in the Fourier domain. The numerical comparison with MCX [18] and experimental comparisons with OCT imaging in phantoms with known optical properties allowed us to assess and validate the performance of MCOCT.

In OCT systems, the high lateral resolution is achieved through tight focusing of the probing beam, which can be performed by the implementation of a Gaussian focusing profile. A method considering the wave nature of photons was developed by focusing the light along the curved trajectories in accordance with Gaussian optics to a diffraction-limited focal spot. In the simulation, a low scattering coefficient was selected to better visualize the curved Gaussian trajectory of photons. When compared to MCX hyperboloid [18], the diffraction-limited spot size at the focus $w_{0,n_2}$ matched well. Optical fluence rates from the two simulators showed an overlap in progression around the beam waist. With distance to the beam waist, the fluence rates varied, which may be due to the different approaches used for mimicking the Gaussian photon trajectory (hyperboloid vs. Gaussian optics). The fluence generated by MCOCT was similarly higher at the surface of the medium ($z=0$ mm) and in-depth ($z=3$ mm) compared to MCX. Using the importance sampling method in MCOCT reduced the number of photons that scattered in the forward direction, which also includes photons that should have exited the bottom boundary of the medium but were biased back to the surface. While this number is reduced for MCOCT, the fluence was corrected by the biased likelihood ratio and therefore does not explain the higher fluences observed with distance from the beam waist. If the likelihood ratios were not corrected for the biased scattering, the disparity between fluences calculated by MCOCT and MCX at the sample’s surface would be different from that of the bottom. When focusing the light in MCOCT, the photon packet propagated by minimal steps allowing us to track the voxel in which the photon packet was located more accurately and adapt the direction in case of changing optical properties. This strategy allowed us to bypass the need to solve an equation of motion previously proposed [30,38]. Then, only $s_\text {min}$ needs to be adapted and mainly depends on the voxel geometry. In MCOCT, the importance sampling method was used to increase the occurrence of backscattering, thereby reducing computational costs for OCT signal simulation, while in MCX all photons scattered according to the HG function. The difference in computation time between the two tools is mainly due to the computational architecture as MCX uses GPU resources while MCOCT uses CPU resources. The performance of the Gaussian illumination was also clear when comparing it with simulations performed with a pencil beam and a geometric focused beam.

When compared with OCT A-lines acquired in phantoms using a short focal length (30 mm), A-lines reconstructed by MCOCT differed within the first $\sim$0.2 mm. This pattern was not observed at a longer focus (75 mm), where MCOCT matched OCT imaging in the four phantoms. This discrepancy is likely due to the reflection of light at the phantom surface, which is not considered in the simulator. When the surface has slight tilts, the coupling of reflected light with the incident light is higher for shorter focal lengths, leading to enhanced distortions in the A-line at the surface. It is also possible that a more accurate detection scheme [30] and scattering function [26,29] would provide a better match between MCOCT and experimental A-lines. The decrease in OCT signals with increasing penetration depth is expected as the probability of photon absorption also increases with greater propagation. A gradual deviation between A-lines generated by MCOCT and phantoms imaging occurred beyond a depth of about 1 to 1.5 mm. This difference is likely due to the fact that fewer photons reached these depths and were detected. The deviation was more prominent in one-layer phantoms imaged at the focal length of 30 mm compared to 75 mm. The reason behind this is unclear to us at present and is a current limitation of the simulator. For OCT applications in biological tissue, this limitation may not be a significant concern as typical OCT systems are designed to provide an imaging depth of approximately $1.5$ mm [6].

The accurate characterization of signal intensity at interfaces is challenging when modeling complex geometry with multiple layers. When the change in refractive index occurs at the interface of two layers, photons undergo either reflection or refraction, which are mathematically described by Fresnel’s equation [6]. MCOCT signal intensity simulated proximal to the interface of the two layers in the bi-layer phantom did not match the one from OCT imaging at either focal lengths. Fresnel’s reflection was not implemented in MCOCT and this discrepancy may be mitigated by considering such reflections as in previous works [30].

Minimizing computational costs is challenging when considering complex media with optical properties mimicking biological tissue. The computation time was higher for phantoms with higher scattering coefficients as these simulations presented more scattering events. For example, generating the OCT signal with the same quantity of detected photons for phantom #1 required higher computational resources than for phantom #3 as scattering was $\sim$6 times higher (see Table 4). Also, the quantity of photon packets launched when initiating the simulation dictates the computation time as well as the ratio of detected photons (0.0525 to 0.33%) that contributed to the OCT signal. As expected, simulations that used $10^9$ photon packets (one-layer phantoms at $f$ = 30 mm) were $\sim$10 times longer than simulations that launched $10^8$ photon packets (one-layer phantoms at $f$ = 75 mm). For the bi-layer phantom, a similar pattern was observed as the simulation using $4\times 10^8$ photon packets (at $f$ = 75 mm) took $\sim$4 times longer than one using $10^8$ photon packets (at $f$ = 30 mm).

In addition to the computational efficiency of the MC process, the number of voxels included in the geometry also affected the computational time of MCOCT. Each time a photon packet propagates by the distance prescribed by $s_\text {min}$, MCOCT estimates its distance to the voxel boundary, which represents several occurrences when considering a geometry with a high number of voxels. This step is essential for geometry designed with multiple layers as optical properties have to be modified when propagating through different tissues. The voxel size also dictates the spatial resolution of the optical fluence rates. A strategy to minimize the number of crossing-boundary events is to use one voxel for each layer with the voxel size in the axial direction being equal to the layer thickness. This simple design can be used for simpler geometry such as rectangular slabs but is less appropriate for complex media. The axial resolution of the simulated OCT data is achieved through the imaging spectrum, whereas the quantity of simulated source/detector positions dictates the lateral resolution of the B-scan. The use of large voxels for simple geometry did not affect the lateral resolution of the OCT image. Another factor that influenced the computation time was the minimal step used to generate the Gaussian focusing beam. The selection of a short minimal step will increase the accuracy of individual photon trajectories but also increase the occurrences with which the trajectory has to be updated. Previous methods have been proposed to automatically adjust the step size during the simulation [38]. In MCOCT, an adjustable step could be implemented but it is not clear how it would reduce computation time while maintaining high accuracy. Another time-consuming procedure in the Monte Carlo module was the random number generation that is needed for each scattering event. Generating one random number each time could be substituted with an a priori array of random numbers generated prior to the simulation.

Our work has several limitations. MCOCT simulations were run on a single CPU thread, which limited speed compared to GPU-based simulators [18,26,28,29]. Accelerating methods such as parallel computing by Compute Unified Device Architecture (CUDA) will be explored for future versions of MCOCT. Using this architecture would also allow performing the initialization, simulation running, postprocessing, and visualization entirely in Matlab. Therefore, users would not need a separate UNIX-emulating console shell to run a compiler and the MC executable, nor would there be a need for intermediate input/output files, which additionally would save storage space. Parallelization would allow simulating a higher quantity of photon packets leading to higher accuracy of the OCT signals at higher sample depths. Noise and/or dispersion were not added to MCOCT signals before reconstruction, but could be added to the reconstructed OCT signal in future works. MCOCT used a scheme of 1000 illumination sources and detectors along the $x$-axis to generate a B-scan image. Since the geometry is described using a 3D voxel system, it is also possible to simulate a scan over the $y$-axis to reconstruct 3D volumes. However, a 3D scan would currently result in a significant increase in computation time. The demonstration of the ability to generate 3D volumes would be more relevant with a code based on parallelization using GPU.

Monte Carlo methods are the gold standard for studying photon propagation in turbid media [7,12,46] and have been instrumental in investigating the relationship between tissue optical properties and imaging results for several techniques including Fourier-domain OCT [22,26]. An accurate analytical/numerical modeling of Fourier-domain OCT imaging is critical to study light-tissue interactions, including speckle contrast [47,48], dispersion [49,50], and multiple scattering [51,52]. A better understanding of these optical phenomena via accurate simulations may help to design an instrument, optimize data acquisition specs and image reconstruction parameters. Numerical and in vitro experiments are often used to produce preliminary results that can guide the development of an in vivo experimental protocol to answer a scientific question.

While further enhancements remain necessary, our open-source simulator can be employed for model validation and system optimization by generating OCT signals across a wavelength spectrum to ascertain the appropriate NA to achieve the desired contrast and its observable depth. Simulating the impact of varying NA or different modes of a photonic lantern [53] can determine the feasibility of achieving a substantial contrast. Our tool can also contribute to the planning of in vivo experiments to provide preliminary estimates of specific biomarkers such as changes in retinal layer thickness to characterize particular physiological/pathological conditions [54]. Previous applications [10,19–25,27,28,30,31] and more recent ones involve the generation of OCT image datasets that can be used in the training of neural networks aimed at tackling the inverse problem in OCT. Some groups already proposed Monte Carlo simulation platforms to generate realistic OCT images for predictive model training [32,33].

5. Conclusion

In this work, a novel Monte Carlo simulator for Fourier-domain OCT (MCOCT) was introduced. The incorporation of a Gaussian light illumination scheme and the ability to increase backscattered event collection allowed us to achieve a more realistic approach to OCT imaging. Optical fluence rates generated by MCOCT were compared with an established MC simulator (MCX), which provided a numerical validation. OCT signals generated by MCOCT were also compared to OCT signals acquired at specific focal lengths in tissue-mimicking phantoms with known optical properties, which provided experimental validation. Future work involves leveraging GPU acceleration for enhanced computational efficiency and speed, as well as considering Fresnel’s reflections proximal to interfaces in multi-layered geometries. MCOCT may further allow the designing of OCT system hardware and image processing for various medical and biomedical applications.

Appendix

Table 5. Variables used to describe Fig. 1 and Fig. 3

View Table | View all tables in this article

 

Fig. 8. Optical coherence tomography (OCT) A-line measured experimentally in Phantom #1 (red) compared to OCT A-lines simulated with a Gaussian (blue), geometric focused (orange), and pencil (green) beam using an imaging lens with a focal length (a) $f=30$ and (b) $f=45$ mm.

Download Full Size | PDF

Funding

Deutscher Akademischer Austauschdienst; Natural Sciences and Engineering Research Council of Canada (#RGPIN-2015-04672, #RGPIN-2018-06151, #RTI-2016-00636).

Acknowledgements

Portions of this work were presented at the SPIE BIOS conference in 2021,116400G in an oral presentation and as a poster at the Optica Biophotonics Congress: Optics in the Life Science in 2023,JTu4B–29 [55,56].

Disclosures

The authors declare no conflicts of interest.

Data Availability

The source code for MCOCT will be available on GitHub upon publication of the manuscript. Raw data used in this work are available upon appropriate request.

References

1. D. Huang, E. A. Swanson, C. P. Lin, et al., “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef]  

2. M. Wojtkowski, R. Leitgeb, A. Kowalczyk, et al., “In vivo human retinal imaging by Fourier domain optical coherence tomography,” J. Biomed. Opt. 7(3), 457–463 (2002). [CrossRef]  

3. I.-K. Jang, B. E. Bouma, D.-H. Kang, et al., “Visualization of coronary atherosclerotic plaques in patients using optical coherence tomography: comparison with intravascular ultrasound,” J. Am. Coll. Cardiol. 39(4), 604–609 (2002). [CrossRef]  

4. T. Gambichler, G. Moussa, M. Sand, et al., “Applications of optical coherence tomography in dermatology,” J. Dermatol. Sci. 40(2), 85–94 (2005). [CrossRef]  

5. F. Benboujja, J. A. Garcia, K. Beaudette, et al., “Intraoperative imaging of pediatric vocal fold lesions using optical coherence tomography,” J. Biomed. Opt. 21(1), 016007 (2016). [CrossRef]  

6. L. V. Wang and H.-I. Wu, Biomedical optics: principles and imaging (John Wiley & Sons, 2012).

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

8. S. T. Flock, M. S. Patterson, B. C. Wilson, et al., “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]  

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

10. 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 (1999). [CrossRef]  

11. D. A. Boas, J. P. Culver, J. J. Stott, et al., “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10(3), 159–170 (2002). [CrossRef]  

12. C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt. 18(5), 050902 (2013). [CrossRef]  

13. S. Jacques, T. Li, and S. Prahl, “mcxyz. c, a 3D Monte Carlo simulation of heterogeneous tissues,” omlc.org/software/mc/mcxyz (2013).

14. D. Marti, R. N. Aasbjerg, P. E. Andersen, et al., “MCmatlab: an open-source, user-friendly, MATLAB-integrated three-dimensional Monte Carlo light transport solver with heat diffusion and tissue damage,” J. Biomed. Opt. 23(12), 1 (2018). [CrossRef]  

15. H. Shen and G. Wang, “A tetrahedron-based inhomogeneous Monte Carlo optical simulator,” Phys. Med. Biol. 55(4), 947–962 (2010). [CrossRef]  

16. Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express 1(1), 165–175 (2010). [CrossRef]  

17. E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt. 13(6), 060504 (2008). [CrossRef]  

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

19. Y. Pan, R. Birngruber, J. Rosperich, et al., “Low-coherence optical tomography in turbid tissue: theoretical analysis,” Appl. Opt. 34(28), 6564–6574 (1995). [CrossRef]  

20. D. J. Smithies, T. Lindmo, Z. Chen, et al., “Signal attenuation and localization in optical coherence tomography studied by Monte Carlo simulation,” Phys. Med. Biol. 43(10), 3025–3044 (1998). [CrossRef]  

21. A. Tycho, T. M. Jørgensen, H. T. Yura, et al., “Derivation of a Monte Carlo method for modeling heterodyne detection in optical coherence tomography systems,” Appl. Opt. 41(31), 6676–6691 (2002). [CrossRef]  

22. Q. Lu, X. Gan, M. Gu, et al., “Monte Carlo modeling of optical coherence tomography imaging through turbid media,” Appl. Opt. 43(8), 1628–1637 (2004). [CrossRef]  

23. G. Xiong, P. Xue, J. Wu, et al., “Particle-fixed Monte Carlo model for optical coherence tomography,” Opt. Express 13(6), 2182–2195 (2005). [CrossRef]  

24. I. T. Lima, A. Kalra, H. E. Hernández-Figueroa, et al., “Fast calculation of multipath diffusive reflectance in optical coherence tomography,” Biomed. Opt. Express 3(4), 692–700 (2012). [CrossRef]  

25. M. Kirillin, I. Meglinski, V. Kuzmin, et al., “Simulation of optical coherence tomography images by Monte Carlo modeling based on polarization vector approach,” Opt. Express 18(21), 21714–21724 (2010). [CrossRef]  

26. J. Mao, Y. Ling, P. Xue, et al., “Monte carlo-based full-wavelength simulator of fourier-domain optical coherence tomography,” Biomed. Opt. Express 13(12), 6317–6334 (2022). [CrossRef]  

27. V. Periyasamy and M. Pramanik, “Importance sampling-based Monte Carlo simulation of time-domain optical coherence tomography with embedded objects,” Appl. Opt. 55(11), 2921–2929 (2016). [CrossRef]  

28. S. Malektaji, I. T. Lima Jr., S. S. Sherif, et al., “Massively parallel simulator of optical coherence tomography of inhomogeneous turbid media,” Comput. Methods Programs Biomed. 150, 97–105 (2017). [CrossRef]  

29. J. Mao, Y. Ling, P. Xue, et al., “Importance sampling-accelerated simulation of full-spectrum backscattered diffuse reflectance,” Biomed. Opt. Express 14(9), 4644–4659 (2023). [CrossRef]  

30. Y. Wang and L. Bai, “Accurate Monte Carlo simulation of frequency-domain optical coherence tomography,” Int. J. for Numer. Methods Biomed. Eng. 35, e3177 (2019). [CrossRef]  

31. A. E. Hartinger, A. S. Nam, I. Chico-Calero, et al., “Monte Carlo modeling of angiographic optical coherence tomography,” Biomed. Opt. Express 5(12), 4338–4349 (2014). [CrossRef]  

32. N. Mavaddat, S. Ahderom, V. Tiporlini, et al., “Simulation of biomedical signals and images using Monte Carlo methods for training of deep learning networks,” in Deep Learning Techniques for Biomedical and Health Informatics, (Elsevier, 2020), pp. 213–236.

33. S. Guo and J. U. Kang, “Convolutional neural network-based common-path optical coherence tomography A-scan boundary-tracking training and validation using a parallel Monte Carlo synthetic dataset,” Opt. Express 30(14), 25876–25890 (2022). [CrossRef]  

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

35. L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]  

36. T. Binzoni, T. S. Leung, A. H. Gandjbakhche, et al., “The use of the Henyey–Greenstein phase function in Monte Carlo simulations in biomedical optics,” Phys. Med. Biol. 51(17), N313–N322 (2006). [CrossRef]  

37. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef]  

38. B. H. Hokr, J. N. Bixler, G. Elpers, et al., “Modeling focusing Gaussian beams in a turbid medium with Monte Carlo simulations,” Opt. Express 23(7), 8699–8705 (2015). [CrossRef]  

39. F. Schneider, J. Draheim, R. Kamberger, et al., “Process and material properties of polydimethylsiloxane (PDMS) for Optical MEMS,” Sensors and Actuators A: Phys. 151(2), 95–99 (2009). [CrossRef]  

40. S. A. Prahl, M. J. van Gemert, and A. J. Welch, “Determining the optical properties of turbid media by using the adding–doubling method,” Appl. Opt. 32(4), 559–568 (1993). [CrossRef]  

41. R. B. Saager, C. Kondru, K. Au, et al., “Multilayer silicone phantoms for the evaluation of quantitative optical techniques in skin imaging,” Proc. SPIE 7567, 756706 (2010). [CrossRef]  

42. J. Schmitt, A. Knüttel, and R. Bonner, “Measurement of optical properties of biological tissues by low-coherence reflectometry,” Appl. Opt. 32(30), 6032–6042 (1993). [CrossRef]  

43. F. J. van der Meer, D. J. Faber, D. M. B. Sassoon, et al., “Localized measurement of optical attenuation coefficients of atherosclerotic plaque constituents by quantitative optical coherence tomography,” IEEE Trans. Med. Imaging 24(10), 1369–1376 (2005). [CrossRef]  

44. M. Almasian, N. Bosschaart, T. G. van Leeuwen, et al., “Validation of quantitative attenuation and backscattering coefficient measurements by optical coherence tomography in the concentration-dependent and multiple scattering regime,” J. Biomed. Opt. 20(12), 121314 (2015). [CrossRef]  

45. A. N. Bashkatov, E. A. Genina, and V. V. Tuchin, “Optical properties of skin, subcutaneous, and muscle tissues: a review,” J. Innov. Opt. Health Sci. 04(01), 9–38 (2011). [CrossRef]  

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

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

48. M. Y. Kirillin, G. Farhat, E. A. Sergeeva, et al., “Speckle statistics in oct images: Monte carlo simulations and experimental studies,” Opt. Lett. 39(12), 3472–3475 (2014). [CrossRef]  

49. V. V. Tuchin, “Tissue optics and photonics: light-tissue interaction,” J. Biomed. Photonics & Eng. 1, 98–134 (2015). [CrossRef]  

50. C. Photiou, E. Bousi, I. Zouvani, et al., “Using speckle to measure tissue dispersion in optical coherence tomography,” Biomed. Opt. Express 8(5), 2528–2535 (2017). [CrossRef]  

51. M. Rinehart, Y. Zhu, and A. Wax, “Quantitative phase spectroscopy,” Biomed. Opt. Express 3(5), 958–965 (2012). [CrossRef]  

52. Y. Ling, Y. Gan, X. Yao, et al., “Phase-noise analysis of swept-source optical coherence tomography systems,” Opt. Lett. 42(7), 1333–1336 (2017). [CrossRef]  

53. M. P. de Sivry-Houle, S. B. Beaudoin, S. Brais-Brunet, et al., “All-fiber few-mode optical coherence tomography using a modally-specific photonic lantern,” Biomed. Opt. Express 12(9), 5704–5719 (2021). [CrossRef]  

54. S. Brais-Brunet, É. Heckel, U. Kanniyappan, et al., “Morphometric and microstructural changes during murine retinal development characterized using in vivo optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 62(13), 20 (2021). [CrossRef]  

55. R. Maltais-Tariant, K. Erdenedalai, U. Kanniyappan, et al., “An open-source platform for optical coherence tomography Monte Carlo simulation,” Proc. SPIE 11640, 116400G (2021). [CrossRef]  

56. K. Erdenedalai, R. Maltais-Taraint, M. Dehaes, et al., “Experimentally Validated Monte Carlo Simulator for Fourier-Domain Optical Coherence Tomography,” in Optical Manipulation and Its Applications, (Optica Publishing Group, 2023), pp. JTu4B–29.