1. Introduction

Optical coherence tomography (OCT) is a non-invasive sub-surface imaging technique that has experienced significant growth in biomedical applications [1, 2]. OCT systems could be implemented with a low-coherence light source and a mechanical scanning sub-system (time-domain OCT). More advanced systems use a low-coherence light source with a spectrometer or a wavelength-swept source (frequency-domain OCT). OCT has an imaging depth that could reach up to 3 mm, depending on the optical properties of the tissue, and it also has one to two orders of magnitude higher resolution than ultrasound imaging. OCT could also produce images inside the body when it is implemented using optical fiber probes. Unlike X-ray or gamma-ray imaging, OCT is safe for biological tissues because it utilizes non-ionizing radiation mainly in the infrared spectrum.

#### 1.1 OCT signal simulation using a Monte Carlo method

The signal obtained by an OCT imaging system consists of ballistic and quasi-ballistic (Class I diffuse reflectance) photons, in addition to multiply scattered photons (Class II diffuse reflectance), that are reflected from tissue [3].

However, multiply scattered photons do not carry practically useful information about the imaged tissue; therefore, they result in a degradation of the OCT signal [4]. In addition, it has been shown that Class II diffuse reflectance represents a fundamental limit related to the imaging depth of OCT in tissue [5]. Therefore, it is important to account for both Class I and Class II diffusive reflectance in any practical simulation of OCT signals.

2. Standard MC method for simulating OCT signals

DOI: http://dx.doi.org/10.5772/intechopen.89555

the Henyey-Greenstein probability density function that is given by

f HGð Þ¼ cosð Þ θ<sup>S</sup>

new propagation direction is θS. Therefore, cos(θS) = u^ �u^<sup>0</sup>

and simulated.

is defined as

5

Our implementation of the MC method to simulate OCT signals is based on Monte Carlo simulation of light transport in multilayered tissues (MCML) [6]. MCML simulates an ensemble of photon packets that are launched as a steady-state pencil beam, normal to the top surface of the medium. Within the tissue, each such photon packet undergoes a random walk whose step size is determined by an exponentially distributed random variable that is parameterized by the interaction coefficient of this tissue. This interaction coefficient is equal to the sum of the absorption μ<sup>a</sup> and the scattering μ<sup>s</sup> coefficients of this tissue. The scattering events that take place at the end of the random steps are characterized by two random angles that determine the next direction of the photon packet. To account for the photon packet scattering, given an anisotropy factor, g, of the tissue, MCML uses

Monte Carlo Methods for Simulation of Optical Coherence Tomography of Turbid Media

<sup>1</sup> � <sup>g</sup><sup>2</sup> 2 1 þ g ð Þ <sup>2</sup> � 2g � cosð Þ θ<sup>S</sup>

where θ<sup>S</sup> is the angle between the propagation direction of the photon packet u^ before the current scattering and u^<sup>0</sup> is the direction of the photon packet after the current scattering. The angle between the previous propagation direction and the

propagation direction is statistically correct, the provisional scattering direction u^<sup>0</sup> is rotated around u^ by an angle ϕ, which is randomly selected from a uniform probability density function with a range from 0 to 2π to generate the propagation direction u^<sup>0</sup> after the current scattering. At the end of each scattering event, the photon packet weight W is reduced according to the step size and the local absorption coefficient μa. The weight W, which is initially set at 1, is proportional to the number of photons in the photon packet. The photon packet is either removed with probability 1/m or is allowed to continue propagating with probability 1�1/m and a new weight equal to <sup>m</sup>�<sup>W</sup> once the weight reaches <sup>W</sup>th <sup>¼</sup> <sup>10</sup>�4. We use the value m = 10 in this work. This procedure, denoted Russian roulette, is an unbiased technique to end simulation of photon packets that have a negligibly low contribution to the Monte Carlo simulation, so that a new photon packet can be initiated

The Class I diffuse reflectance at depth equal to z is obtained by calculating the mean value of an indicator function I<sup>1</sup> that represents a spatial and temporal

indicator function I<sup>1</sup> of such spatial and temporal filter for the ith photon packet

<sup>I</sup>1ð Þ¼ <sup>z</sup>, <sup>i</sup> 1, ri <sup>&</sup>lt; <sup>d</sup>max, <sup>θ</sup><sup>z</sup>,<sup>i</sup> <sup>&</sup>lt; <sup>θ</sup>max, j j <sup>Δ</sup>Si � <sup>2</sup><sup>z</sup> <sup>&</sup>lt;lc

where lc is the optical source's coherence length; ri is the distance to the origin of the ith reflected photon packet; dmax and θmax are the maximum photon packet collecting diameter and angle, respectively; θ<sup>z</sup>,<sup>i</sup> is the angle with the z-axis (normal to the tissue interface); Δsi is the optical path; and z is the photon packet's maximum depth.

At any depth, the diffuse reflectance R<sup>1</sup> is the expected value of I<sup>1</sup> at this depth,

filter of Class I diffuse reflectance for all simulated photon packets. The

0, othwerise �

R1ð Þ¼ z

1 N X N

i¼1

and its standard deviation σR,1 could be estimated by

<sup>3</sup>=<sup>2</sup> , (1)

. To ensure that the new

, (2)

I1ð Þ z, i W ið Þ (3)

Since it is not practical to simulate light transport in turbid media, for example, tissue, using electromagnetic waves, especially due to diffusive scattering, a Monte Carlo (MC) method of simulating light transport in tissue has been typically used [6–9]. However, the computational cost of this MC-based simulation of OCT systems could be very high, as the probability of detecting diffusively reflected photons from tissue is very low [4, 5].

To reduce the computational cost, thereby accelerate, this MC simulation, importance sampling could be used to speed up simulations by orders of magnitude. Importance sampling has been applied earlier to optical communications [10, 11], confocal microscopy [12], atmospheric optics [13], and diffuse optical tomography [8].

To improve the computational efficiency of the MC-based simulation of OCT systems [6], Yao and Wang proposed the first importance sampling technique to simulate OCT signals from a multilayered turbid medium [3]. However, their method only enabled the simulation of OCT signals from a thin shallow layer in tissue, as the results obtained from deeper tissue regions were underestimated. In [14], we, the authors of this chapter, developed another more advanced importance sampling technique by implementing multiple biased scatterings per photon packet, and by developing a photon splitting procedure. Our advanced importance sampling resulted in more accurate and computationally efficient simulations of OCT signals due to ballistic and quasi-ballistic photons. However, it still underestimated OCT signals due to multiply scattered photons. To enable accurate simulation of OCT signals due to both Class I and Class II diffusive reflectances, we further developed our importance sampling technique to accurately and efficiently simulate diffusive reflectance due to photons that undergo multiple scattering [15]. In this method, additional biased scatterings were randomly applied, which enabled accurate simulation of both Class I and Class II diffusive reflectances, with a speed-up of three orders of magnitude compared to the standard MC method.

Our advanced importance sampling techniques above were implemented to simulate OCT of tissues with planar geometry [6]. To enable simulation of OCT of arbitrarily shaped turbid media, we used the mesh-based MC method of light transport in tissue proposed by Fang [16]. This method uses a Plücker coordinate system to efficiently calculate intersections between paths of light propagation with interfaces of the object regions that are modeled using tetrahedrons. We combined this mesh-based MC method with our importance sampling techniques to simulate OCT signals from tissue with arbitrarily shaped regions. However, since it was still computationally costly to simulate a full OCT B-scan using this method, we also developed a parallel implementation of this simulator that exploited the massively parallel computing capabilities of Graphics Processing Units (GPUs) to accelerate this simulator by two additional orders of magnitude [17, 18]. This GPU-based implementation enabled simulation of OCT B-scans of arbitrarily shaped turbid media in a few minutes using commonly available workstations.

In Section 2 of this chapter, we present a standard MC method for simulating OCT signals. In Section 3, we present our first importance sampling implementation that enables the simulation of OCT signals from higher depths inside turbid media. In Section 4, we present our more advanced importance sampling implementation that accurately calculates both Class I and Class II diffusive reflectances, and is three orders of magnitude faster than the standard MC simulation.

Monte Carlo Methods for Simulation of Optical Coherence Tomography of Turbid Media DOI: http://dx.doi.org/10.5772/intechopen.89555

## 2. Standard MC method for simulating OCT signals

However, multiply scattered photons do not carry practically useful information about the imaged tissue; therefore, they result in a degradation of the OCT signal [4]. In addition, it has been shown that Class II diffuse reflectance represents a fundamental limit related to the imaging depth of OCT in tissue [5]. Therefore, it is important to account for both Class I and Class II diffusive reflectance in any

Theory, Application, and Implementation of Monte Carlo Method in Science and Technology

Since it is not practical to simulate light transport in turbid media, for example, tissue, using electromagnetic waves, especially due to diffusive scattering, a Monte Carlo (MC) method of simulating light transport in tissue has been typically used [6–9]. However, the computational cost of this MC-based simulation of OCT systems could be very high, as the probability of detecting diffusively reflected pho-

To reduce the computational cost, thereby accelerate, this MC simulation, importance sampling could be used to speed up simulations by orders of magnitude. Importance sampling has been applied earlier to optical communications [10, 11],

To improve the computational efficiency of the MC-based simulation of OCT systems [6], Yao and Wang proposed the first importance sampling technique to simulate OCT signals from a multilayered turbid medium [3]. However, their method only enabled the simulation of OCT signals from a thin shallow layer in tissue, as the results obtained from deeper tissue regions were underestimated. In [14], we, the authors of this chapter, developed another more advanced importance sampling technique by implementing multiple biased scatterings per photon packet, and by developing a photon splitting procedure. Our advanced importance sampling resulted in more accurate and computationally efficient simulations of OCT signals due to ballistic and quasi-ballistic photons. However, it still underestimated OCT signals due to multiply scattered photons. To enable accurate simulation of OCT signals due to both Class I and Class II diffusive reflectances, we further developed our importance sampling technique to accurately and efficiently simulate diffusive reflectance due to photons that undergo multiple scattering [15]. In this method, additional biased scatterings were randomly applied, which enabled accurate simulation of both Class I and Class II diffusive reflectances, with a speed-up of

confocal microscopy [12], atmospheric optics [13], and diffuse optical

three orders of magnitude compared to the standard MC method.

media in a few minutes using commonly available workstations.

orders of magnitude faster than the standard MC simulation.

Our advanced importance sampling techniques above were implemented to simulate OCT of tissues with planar geometry [6]. To enable simulation of OCT of arbitrarily shaped turbid media, we used the mesh-based MC method of light transport in tissue proposed by Fang [16]. This method uses a Plücker coordinate system to efficiently calculate intersections between paths of light propagation with interfaces of the object regions that are modeled using tetrahedrons. We combined this mesh-based MC method with our importance sampling techniques to simulate OCT signals from tissue with arbitrarily shaped regions. However, since it was still computationally costly to simulate a full OCT B-scan using this method, we also developed a parallel implementation of this simulator that exploited the massively parallel computing capabilities of Graphics Processing Units (GPUs) to accelerate this simulator by two additional orders of magnitude [17, 18]. This GPU-based implementation enabled simulation of OCT B-scans of arbitrarily shaped turbid

In Section 2 of this chapter, we present a standard MC method for simulating OCT signals. In Section 3, we present our first importance sampling implementation that enables the simulation of OCT signals from higher depths inside turbid media. In Section 4, we present our more advanced importance sampling implementation that accurately calculates both Class I and Class II diffusive reflectances, and is three

practical simulation of OCT signals.

tons from tissue is very low [4, 5].

tomography [8].

4

Our implementation of the MC method to simulate OCT signals is based on Monte Carlo simulation of light transport in multilayered tissues (MCML) [6]. MCML simulates an ensemble of photon packets that are launched as a steady-state pencil beam, normal to the top surface of the medium. Within the tissue, each such photon packet undergoes a random walk whose step size is determined by an exponentially distributed random variable that is parameterized by the interaction coefficient of this tissue. This interaction coefficient is equal to the sum of the absorption μ<sup>a</sup> and the scattering μ<sup>s</sup> coefficients of this tissue. The scattering events that take place at the end of the random steps are characterized by two random angles that determine the next direction of the photon packet. To account for the photon packet scattering, given an anisotropy factor, g, of the tissue, MCML uses the Henyey-Greenstein probability density function that is given by

$$f\_{HG}(\cos(\theta\_{\mathcal{S}})) = \frac{\mathbf{1} - \mathbf{g}^2}{2(\mathbf{1} + \mathbf{g}^2 - 2\mathbf{g} \cdot \cos(\theta\_{\mathcal{S}}))^{3/2}},\tag{1}$$

where θ<sup>S</sup> is the angle between the propagation direction of the photon packet u^ before the current scattering and u^<sup>0</sup> is the direction of the photon packet after the current scattering. The angle between the previous propagation direction and the new propagation direction is θS. Therefore, cos(θS) = u^ �u^<sup>0</sup> . To ensure that the new propagation direction is statistically correct, the provisional scattering direction u^<sup>0</sup> is rotated around u^ by an angle ϕ, which is randomly selected from a uniform probability density function with a range from 0 to 2π to generate the propagation direction u^<sup>0</sup> after the current scattering. At the end of each scattering event, the photon packet weight W is reduced according to the step size and the local absorption coefficient μa. The weight W, which is initially set at 1, is proportional to the number of photons in the photon packet. The photon packet is either removed with probability 1/m or is allowed to continue propagating with probability 1�1/m and a new weight equal to <sup>m</sup>�<sup>W</sup> once the weight reaches <sup>W</sup>th <sup>¼</sup> <sup>10</sup>�4. We use the value m = 10 in this work. This procedure, denoted Russian roulette, is an unbiased technique to end simulation of photon packets that have a negligibly low contribution to the Monte Carlo simulation, so that a new photon packet can be initiated and simulated.

The Class I diffuse reflectance at depth equal to z is obtained by calculating the mean value of an indicator function I<sup>1</sup> that represents a spatial and temporal filter of Class I diffuse reflectance for all simulated photon packets. The indicator function I<sup>1</sup> of such spatial and temporal filter for the ith photon packet is defined as

$$I\_1(\mathbf{z}, i) = \begin{cases} \mathbf{1}, & r\_i < d\_{\text{max}}, \quad \theta\_{\mathbf{z}, i} < \theta\_{\text{max}}, \quad |\Delta \mathbf{S}\_i - \mathbf{2}\mathbf{z}| < l\_c \\\ \mathbf{0}, & \text{otherwise} \end{cases},\tag{2}$$

where lc is the optical source's coherence length; ri is the distance to the origin of the ith reflected photon packet; dmax and θmax are the maximum photon packet collecting diameter and angle, respectively; θ<sup>z</sup>,<sup>i</sup> is the angle with the z-axis (normal to the tissue interface); Δsi is the optical path; and z is the photon packet's maximum depth.

At any depth, the diffuse reflectance R<sup>1</sup> is the expected value of I<sup>1</sup> at this depth, and its standard deviation σR,1 could be estimated by

$$R\_1(\mathbf{z}) = \frac{1}{N} \sum\_{i=1}^{N} I\_1(\mathbf{z}, i) W(i) \tag{3}$$

and

$$\sigma\_{R,1}^2(\mathbf{z}) = \frac{1}{N(N-1)} \sum\_{i=1}^{N} \left[ I\_1(\mathbf{z}, i) W(i) - R\_1 \right]^2 \tag{4}$$

uniform probability density function with a range from 0 to 2π to generate the propagation direction u^<sup>0</sup> after the current scattering. This procedure ensures a more accurate model of the light scattering in tissue. The difference in the rotation by ϕ between the model with importance sampling and the standard model is that the rotation in the model with importance sampling is done around the biased direction v^, while the rotation in the standard model in is done around the direction u^ prior to the current scattering. After the first biased scattering, this procedure produces the new propagation direction u^<sup>0</sup> of the photon packet. Afterward, the scattered photon packet is associated with a likelihood ratio as discussed in other applications of this method [10, 11]. Using our biased angle's probability density function, the likeli-

Monte Carlo Methods for Simulation of Optical Coherence Tomography of Turbid Media

<sup>¼</sup> <sup>1</sup> � <sup>g</sup><sup>2</sup> 1 � a<sup>2</sup>

where cos(θS) = u^ �u^<sup>0</sup> is determined, after the biased scattering, from the randomly sampled values of θ<sup>B</sup> and ϕ. The ratio of the probability of the scattering angle appearing in the biased case with the standard case is the likelihood ratio that is shown in Eq. (6). In addition to depending on θB, the likelihood ratio also depends also on θS. Figure 1 shows a schematic drawing of these vectors and the angles used in this direction biasing procedure. Note that the choice of bias distribution only affects the speed of convergence of the simulation. Therefore, other biased probability function could also be used to randomly generate the biased scattering toward

As a photon packet is biased toward the apparent position of the collecting optical fiber, at any given depth in the tissue, the photon packet becomes more likely to be collected at the tip of the fiber. However, the photon packet could be scattered several times after the first backscatter bias before reaching the optical collection system. These additional scatterings, according to Eq. (1), reduce the correlation between the biased direction and the event in which the photon packet is collected. We overcome this reduction in correlation by continuing to bias the

Schematic representation of vectors and angles used in biasing the scattering direction. Reprinted with

permission from [14] © The Optical Society of America.

<sup>1</sup> <sup>þ</sup> <sup>a</sup><sup>2</sup> � <sup>2</sup><sup>a</sup> � cosð Þ <sup>θ</sup><sup>B</sup> 1 þ g<sup>2</sup> � 2g � cosð Þ θ<sup>S</sup> <sup>3</sup>=<sup>2</sup>

, (6)

hood ratio of the photon packet, Eq. (5), is given by

DOI: http://dx.doi.org/10.5772/intechopen.89555

f HGð Þ cosð Þ θ<sup>S</sup> f <sup>B</sup>ð Þ cosð Þ θ<sup>B</sup>

3.2 Scattering angles of additional events of backscattering

Lð Þ¼ cosð Þ θ<sup>B</sup>

the bias direction v^.

Figure 1.

7

where N is the number of photons packets used in MC-based simulations.

### 3. Importance sampling for simulation of Class I OCT signal

Our first importance sampling technique to simulate OCT signals aimed at increasing the number of photons collected at the detector. This algorithm uses the same method described in Section 2, where we also use the same square time gating given by [3].

Since most tissues are highly forward-scattering, their anisotropy factor is close to 1. Therefore, there is a very small probability that a simulated photon packet at any given depth in the tissue would undergo scattering in the backward direction toward the OCT probe. The probability of collecting Class I photons drops rapidly with depth in the tissue from which the photon is scattered in the backward direction. To allow faster simulation of Class I photons, we designed an importance sampling method that biases the direction u^<sup>0</sup> of a scattered photon toward the tip of the light-collecting optical fiber, v^, as the photon packet reaches the depth range of interest. By defining the origin of a Cartesian coordinate system at the center of the tip of the light-collecting optical fiber, the bias direction in which this fiber is located is defined as v^ ¼ �R=j j R , where R ¼ xx^ þ y^y þ zz^ is the position vector of the scattering location in the tissue.

All photon packets propagating in a direction close to v^ will contribute to the simulated Class I diffuse reflectance with a higher probability. Therefore, this bias direction is more efficient than biasing only in the backward direction, which may not be consistent with the direction of the light-collecting optical fiber. This choice of the bias direction is particularly effective for photon packets propagating deep in the tissue, where such photon packets experience one or more scattering events before they are diffusively reflected.

#### 3.1 Scattering angle due to first event of backscattering

As the photon packet reaches the depth range targeted, the propagation direction u^<sup>0</sup> of the scattered photon packet is biased toward the bias direction v^, as opposed to being most likely scattered close to the previous propagation direction u^ as in the practical case with anisotropy g close to 1 and different from the bias toward �u, ^ the opposite of the direction of propagation, as it is done in [3]. To randomly select the biased angle θ<sup>B</sup> between the new scattering direction u^<sup>0</sup> and the biased direction v, ^ we use the same probability density function in Eq. (1). However, the bias coefficient does not have to correspond to the anisotropy factor g. Therefore, the probability density function of the biased angle is given by

$$f\_B(\cos(\theta\_B)) = \frac{1 - a^2}{2(1 + a^2 - 2a \cdot \cos(\theta\_B))^{3/2}},\tag{5}$$

where a is a bias coefficient. After randomly sampling a biased angle θ<sup>B</sup> away from the biased direction v^, so that cos(θB) = v^ �u^<sup>0</sup> , the provisional scattering direction u^<sup>0</sup> is rotated around v^ by an angle ϕ, which is randomly selected from a Monte Carlo Methods for Simulation of Optical Coherence Tomography of Turbid Media DOI: http://dx.doi.org/10.5772/intechopen.89555

uniform probability density function with a range from 0 to 2π to generate the propagation direction u^<sup>0</sup> after the current scattering. This procedure ensures a more accurate model of the light scattering in tissue. The difference in the rotation by ϕ between the model with importance sampling and the standard model is that the rotation in the model with importance sampling is done around the biased direction v^, while the rotation in the standard model in is done around the direction u^ prior to the current scattering. After the first biased scattering, this procedure produces the new propagation direction u^<sup>0</sup> of the photon packet. Afterward, the scattered photon packet is associated with a likelihood ratio as discussed in other applications of this method [10, 11]. Using our biased angle's probability density function, the likelihood ratio of the photon packet, Eq. (5), is given by

$$L(\cos(\theta\_B)) = \frac{f\_{HG}(\cos(\theta\_S))}{f\_B(\cos(\theta\_B))} = \frac{1-\mathbf{g}^2}{\mathbf{1}-a^2} \left(\frac{\mathbf{1}+a^2-2a\cdot\cos(\theta\_B)}{\mathbf{1}+\mathbf{g}^2-2\mathbf{g}\cdot\cos(\theta\_S)}\right)^{3/2},\tag{6}$$

where cos(θS) = u^ �u^<sup>0</sup> is determined, after the biased scattering, from the randomly sampled values of θ<sup>B</sup> and ϕ. The ratio of the probability of the scattering angle appearing in the biased case with the standard case is the likelihood ratio that is shown in Eq. (6). In addition to depending on θB, the likelihood ratio also depends also on θS. Figure 1 shows a schematic drawing of these vectors and the angles used in this direction biasing procedure. Note that the choice of bias distribution only affects the speed of convergence of the simulation. Therefore, other biased probability function could also be used to randomly generate the biased scattering toward the bias direction v^.

#### 3.2 Scattering angles of additional events of backscattering

As a photon packet is biased toward the apparent position of the collecting optical fiber, at any given depth in the tissue, the photon packet becomes more likely to be collected at the tip of the fiber. However, the photon packet could be scattered several times after the first backscatter bias before reaching the optical collection system. These additional scatterings, according to Eq. (1), reduce the correlation between the biased direction and the event in which the photon packet is collected. We overcome this reduction in correlation by continuing to bias the

#### Figure 1.

Schematic representation of vectors and angles used in biasing the scattering direction. Reprinted with permission from [14] © The Optical Society of America.

and

given by [3].

σ2 <sup>R</sup>,1ð Þ¼ z

the scattering location in the tissue.

before they are diffusively reflected.

3.1 Scattering angle due to first event of backscattering

f <sup>B</sup>ð Þ¼ cosð Þ θ<sup>B</sup>

from the biased direction v^, so that cos(θB) = v^ �u^<sup>0</sup>

6

1 N Nð Þ � 1

3. Importance sampling for simulation of Class I OCT signal

X N

½ � I1ð Þ z, i W iðÞ� R<sup>1</sup>

<sup>2</sup> (4)

i¼1

where N is the number of photons packets used in MC-based simulations.

Theory, Application, and Implementation of Monte Carlo Method in Science and Technology

Our first importance sampling technique to simulate OCT signals aimed at increasing the number of photons collected at the detector. This algorithm uses the same method described in Section 2, where we also use the same square time gating

Since most tissues are highly forward-scattering, their anisotropy factor is close to 1. Therefore, there is a very small probability that a simulated photon packet at any given depth in the tissue would undergo scattering in the backward direction toward the OCT probe. The probability of collecting Class I photons drops rapidly with depth in the tissue from which the photon is scattered in the backward direction. To allow faster simulation of Class I photons, we designed an importance sampling method that biases the direction u^<sup>0</sup> of a scattered photon toward the tip of the light-collecting optical fiber, v^, as the photon packet reaches the depth range of interest. By defining the origin of a Cartesian coordinate system at the center of the tip of the light-collecting optical fiber, the bias direction in which this fiber is located is defined as v^ ¼ �R=j j R , where R ¼ xx^ þ y^y þ zz^ is the position vector of

All photon packets propagating in a direction close to v^ will contribute to the simulated Class I diffuse reflectance with a higher probability. Therefore, this bias direction is more efficient than biasing only in the backward direction, which may not be consistent with the direction of the light-collecting optical fiber. This choice of the bias direction is particularly effective for photon packets propagating deep in the tissue, where such photon packets experience one or more scattering events

As the photon packet reaches the depth range targeted, the propagation direction u^<sup>0</sup> of the scattered photon packet is biased toward the bias direction v^, as opposed to being most likely scattered close to the previous propagation direction u^ as in the practical case with anisotropy g close to 1 and different from the bias toward �u, ^ the opposite of the direction of propagation, as it is done in [3]. To randomly select the biased angle θ<sup>B</sup> between the new scattering direction u^<sup>0</sup> and the biased direction v, ^ we use the same probability density function in Eq. (1). However, the bias coefficient does not have to correspond to the anisotropy factor g. Therefore, the probability density function of the biased angle is given by

> <sup>1</sup> � <sup>a</sup><sup>2</sup> 2 1 þ a ð Þ <sup>2</sup> � 2a � cosð Þ θ<sup>B</sup>

where a is a bias coefficient. After randomly sampling a biased angle θ<sup>B</sup> away

direction u^<sup>0</sup> is rotated around v^ by an angle ϕ, which is randomly selected from a

<sup>3</sup>=<sup>2</sup> , (5)

, the provisional scattering

scattering direction u^<sup>0</sup> toward the direction v^, pointing to the apparent position of the optical collection system, at every scattering point until the photon packet is removed. These additional biases still use both Eqs. (7) and (8). Since the random values drawn for the angle between the scattering direction and the biased direction are independent of each other and are also independent of the previous scattering events, the overall likelihood ratio of a collected photon packet results from the multiplication of all the likelihood ratios of all the biased scattering in that particular simulation.

Once a photon packet experiences the first biased scattering, that photon packet is biased at all additional scattering points until it is removed from the simulation, which can occur when the photon packet is removed by Russian roulette, as described in Section 2, or it leaves the tissue. After simulating N launched photon packets using importance sampling, the diffuse reflectance R<sup>1</sup> and its standard deviation σR,1 could be calculated with

$$R\_1(\mathbf{z}) = \frac{1}{N} \sum\_{i=1}^{N} I\_1(\mathbf{z}, i) L(i) W(i) \tag{7}$$

L ið Þ< 1, also undergoes biased backscattering in the tissue at the end of the next step, which could result in another photon packet split, and successive additional biased scatterings toward the tip of the collecting optical fiber until the photon packet propagates beyond the simulation domain. In cases that we investigated, this procedure increased the computational time of each photon packet by five times when compared with a photon packet computed using the standard Monte Carlo method with the same number of launched photon packets N. The increase in the computational time of our importance sampling-based implementation, compared to the standard method with the same number of launched photon packets, depends on the average value of the mean free path, and on the width of the target depth range. We note that in our importance sampling implementation, we do not count split photon packets as additional photon packets when determining the value of N in Eqs. (7) and (8), as the use of their corresponding likelihood ratios will generate the correct result. As a photon packet propagates beyond the target region, the packet will propagate using the standard scattering procedure until it is terminated. Once this photon packet is terminated, a new photon packet will be simulated from the OCT probe, as it is the case in the standard MCML. Even though the splitting procedure implies that the cost of simulating a launched photon with this importance sampling method is higher than the computational cost of simulating launched photos using the standard MCML, the computation cost of the Class I diffuse reflectance in our Monte Carlo simulations with importance sampling required as little as one-thousandth of the computational cost required by the standard Monte Carlo method to achieve the same accuracy in the calculated diffuse

Monte Carlo Methods for Simulation of Optical Coherence Tomography of Turbid Media

DOI: http://dx.doi.org/10.5772/intechopen.89555

We validate our importance sampling technique for simulation of OCT signals from multilayered tissue, with different refractive indices and scattering properties, by comparing its results with those obtained by the standard Monte Carlo method. As shown in Figure 2, light is emitted by an optical fiber probe that is reflected by

The shown optical system has a focusing lens with a numerical aperture (NA) that allows collecting light at an angle of up to 4° and a diameter of 0.5 mm. Similar to the setups in [5, 10, 11], we assume a point source that emits in the vertical direction. Air is present between the center of the probe and the first layer of tissue, which is placed 2.12 mm from the center of the fiber. We simulate a three-layer turbid medium with refractive-index mismatch at its interfaces. The first layer, extending from 2.12 to 2.22 mm from the tip of the fiber, has absorption coefficient

, and refractive index n = 1.

, scattering coefficient μ<sup>s</sup> = 60 cm�<sup>1</sup>

Schematic representation of a setup to simulate OCT signals. Reprinted with permission from [14]

reflectance.

a prism.

μ<sup>a</sup> = 1.5 cm�<sup>1</sup>

Figure 2.

9

© The Optical Society of America.

3.4 Numerical results

and

$$\sigma\_{R,1}^2(z) = \frac{1}{N(N-1)} \sum\_{i=1}^{N} \left[ I\_1(z,i)L(i)\mathcal{W}(i) - R\_1 \right]^2. \tag{8}$$

Eqs. (7) and (8) are similar to Eqs. (3) and (4), except that the indicator function is multiplied by its corresponding likelihood ratio. Using this method, a significantly larger number of photon packets are scattered from a specific depth range toward the collecting optical system than the number obtained using a standard MCML implementation. At the end of this biased simulation, each photon packet is weighted by its likelihood ratio, which adjusts the contribution of each packet to the estimation of the Class I diffuse reflectance. The estimated diffuse reflectance converges toward its true value faster, by several orders of magnitude, when compared to the standard Monte Carlo method.

#### 3.3 Importance sampling effectiveness and depth of tissue

One drawback of previously existing bias methods, for example, [3, 7, 8] is an underestimation of the diffuse reflectance beyond the targeted depth range. The application of the first backward bias reduces the probability that this photon packet would propagate beyond that portion of the tissue. This would lead to a statistical bias to this importance sampling method similar to that in the angle biasing procedure used in [3] and the method used in [7, 8], which limits the effectiveness of those methods to a thin target layer.

We make sure we obtain correct statistics by splitting the photon packet into two photon packets before the first biased scattering [14]. The first of these two photon packets is the one biased toward the collecting optical system. The second photon packet starts propagating from the location in which the biased backscattering occurred, where its initial direction calculated by applying the standard procedure to the previous direction u^ as shown in Section 3.1. To ensure that there is no statistical bias associated with the forward-propagating photon packet that was split, it will be assigned a likelihood ratio L<sup>0</sup> ð Þi , which is a complement to the likelihood of the biased backward scattered photon packet L ið Þ such that L<sup>0</sup> ðÞ¼ i 1 � L ið Þ to this second photon packet. This second photon packet, only generated if Monte Carlo Methods for Simulation of Optical Coherence Tomography of Turbid Media DOI: http://dx.doi.org/10.5772/intechopen.89555

L ið Þ< 1, also undergoes biased backscattering in the tissue at the end of the next step, which could result in another photon packet split, and successive additional biased scatterings toward the tip of the collecting optical fiber until the photon packet propagates beyond the simulation domain. In cases that we investigated, this procedure increased the computational time of each photon packet by five times when compared with a photon packet computed using the standard Monte Carlo method with the same number of launched photon packets N. The increase in the computational time of our importance sampling-based implementation, compared to the standard method with the same number of launched photon packets, depends on the average value of the mean free path, and on the width of the target depth range. We note that in our importance sampling implementation, we do not count split photon packets as additional photon packets when determining the value of N in Eqs. (7) and (8), as the use of their corresponding likelihood ratios will generate the correct result. As a photon packet propagates beyond the target region, the packet will propagate using the standard scattering procedure until it is terminated. Once this photon packet is terminated, a new photon packet will be simulated from the OCT probe, as it is the case in the standard MCML. Even though the splitting procedure implies that the cost of simulating a launched photon with this importance sampling method is higher than the computational cost of simulating launched photos using the standard MCML, the computation cost of the Class I diffuse reflectance in our Monte Carlo simulations with importance sampling required as little as one-thousandth of the computational cost required by the standard Monte Carlo method to achieve the same accuracy in the calculated diffuse reflectance.

#### 3.4 Numerical results

scattering direction u^<sup>0</sup> toward the direction v^, pointing to the apparent position of the optical collection system, at every scattering point until the photon packet is removed. These additional biases still use both Eqs. (7) and (8). Since the random values drawn for the angle between the scattering direction and the biased direction are independent of each other and are also independent of the previous scattering events, the overall likelihood ratio of a collected photon packet results from the multiplication of all the likelihood ratios of all the biased scattering in that particular

Theory, Application, and Implementation of Monte Carlo Method in Science and Technology

Once a photon packet experiences the first biased scattering, that photon packet is biased at all additional scattering points until it is removed from the simulation, which can occur when the photon packet is removed by Russian roulette, as described in Section 2, or it leaves the tissue. After simulating N launched photon packets using importance sampling, the diffuse reflectance R<sup>1</sup> and its standard

I1ð Þ z, i L ið ÞW ið Þ (7)

2

ð Þi , which is a complement to the

ðÞ¼ i

: (8)

½ � I1ð Þ z, i L ið ÞW iðÞ� R<sup>1</sup>

simulation.

and

deviation σR,1 could be calculated with

σ2 <sup>R</sup>,1ð Þ¼ z

R1ð Þ¼ z

1 N Nð Þ � 1

when compared to the standard Monte Carlo method.

effectiveness of those methods to a thin target layer.

split, it will be assigned a likelihood ratio L<sup>0</sup>

8

3.3 Importance sampling effectiveness and depth of tissue

1 N X N

i¼1

X N

i¼1

Eqs. (7) and (8) are similar to Eqs. (3) and (4), except that the indicator function is multiplied by its corresponding likelihood ratio. Using this method, a significantly larger number of photon packets are scattered from a specific depth range toward the collecting optical system than the number obtained using a standard MCML implementation. At the end of this biased simulation, each photon packet is weighted by its likelihood ratio, which adjusts the contribution of each packet to the estimation of the Class I diffuse reflectance. The estimated diffuse reflectance converges toward its true value faster, by several orders of magnitude,

One drawback of previously existing bias methods, for example, [3, 7, 8] is an underestimation of the diffuse reflectance beyond the targeted depth range. The application of the first backward bias reduces the probability that this photon packet would propagate beyond that portion of the tissue. This would lead to a statistical bias to this importance sampling method similar to that in the angle biasing procedure used in [3] and the method used in [7, 8], which limits the

We make sure we obtain correct statistics by splitting the photon packet into two photon packets before the first biased scattering [14]. The first of these two photon packets is the one biased toward the collecting optical system. The second photon packet starts propagating from the location in which the biased backscattering occurred, where its initial direction calculated by applying the standard procedure to the previous direction u^ as shown in Section 3.1. To ensure that there is no statistical bias associated with the forward-propagating photon packet that was

likelihood of the biased backward scattered photon packet L ið Þ such that L<sup>0</sup>

1 � L ið Þ to this second photon packet. This second photon packet, only generated if

We validate our importance sampling technique for simulation of OCT signals from multilayered tissue, with different refractive indices and scattering properties, by comparing its results with those obtained by the standard Monte Carlo method. As shown in Figure 2, light is emitted by an optical fiber probe that is reflected by a prism.

The shown optical system has a focusing lens with a numerical aperture (NA) that allows collecting light at an angle of up to 4° and a diameter of 0.5 mm. Similar to the setups in [5, 10, 11], we assume a point source that emits in the vertical direction. Air is present between the center of the probe and the first layer of tissue, which is placed 2.12 mm from the center of the fiber. We simulate a three-layer turbid medium with refractive-index mismatch at its interfaces. The first layer, extending from 2.12 to 2.22 mm from the tip of the fiber, has absorption coefficient μ<sup>a</sup> = 1.5 cm�<sup>1</sup> , scattering coefficient μ<sup>s</sup> = 60 cm�<sup>1</sup> , and refractive index n = 1.

#### Figure 2.

Schematic representation of a setup to simulate OCT signals. Reprinted with permission from [14] © The Optical Society of America.

f <sup>B</sup>ð Þ¼ cosð Þ θ<sup>B</sup>

given by

u^0 .

formula

11

cos <sup>θ</sup><sup>B</sup>,<sup>i</sup> <sup>¼</sup> <sup>1</sup>

2a

ing is calculated according to the formula

Lð Þ¼ cos θ<sup>B</sup>

<sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>1</sup> � ui

4.2 Scattering angles of additional biased backscatterings

Lð Þ¼ cosð Þ θ<sup>B</sup>

<sup>1</sup> � <sup>1</sup>�<sup>a</sup> ffiffiffiffiffiffiffi <sup>a</sup>2þ<sup>1</sup> <sup>p</sup>

8 ><

DOI: http://dx.doi.org/10.5772/intechopen.89555

>:

OCT probe v^, where cosð Þ¼ θ<sup>B</sup> v^ � u^<sup>0</sup>

f HGð Þ cosð Þ θ<sup>S</sup> f <sup>B</sup>ð Þ cosð Þ θ<sup>B</sup>

0, otherwise

� ��<sup>1</sup> að Þ 1 � a

Monte Carlo Methods for Simulation of Optical Coherence Tomography of Turbid Media

<sup>¼</sup> <sup>1</sup> � <sup>g</sup><sup>2</sup> 2að Þ 1 � a

1 þ a ð Þ <sup>2</sup> � 2a cosð Þ θ<sup>B</sup>

<sup>1</sup> � <sup>1</sup> � <sup>a</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>1</sup> <sup>p</sup>

where cosð Þ¼ θ<sup>S</sup> u^ � u^<sup>0</sup> . We note that cosð Þ θ<sup>B</sup> is obtained using the probability density function in Eq. (9), where it is used to obtain the new propagation direction

To sample angles according to the biased probability density function in (8), one could use any uniform pseudo-random number generator that would be typically available in scientific software libraries. For example, if ui is a random number distributed uniformly between 0 and 1, a random value for cosð Þ θ<sup>B</sup> that satisfies Eq. (9) with bias coefficient a could be generated with the following conversion

> <sup>1</sup> � <sup>a</sup> � <sup>1</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>1</sup> <sup>p</sup>

A second enhancement that could be made to the importance sampling technique, described in Section 3, is to bias the additional scatterings toward the center of the OCT probe v^ with probability 0 ≤ p ≤ 1. That contrasts with the technique presented in Section 3, in which p = 1 (all the additional scatterings were biased). An unbiased scattering is performed in case a bias scattering is not applied in a given point where scattering takes place. The likelihood ratio associated with this scatter-

If the biased function f <sup>B</sup>ð Þ cosð Þ θ<sup>B</sup> is selected to sample a random value of cosð Þ θ<sup>S</sup> , which is an event with probability p, cosð Þ¼ θ<sup>S</sup> u^ � u^<sup>0</sup> is a function of cosð Þ θ<sup>B</sup> that is statistically sampled from the probability density function in Eq. (9).

f HGð Þ cos θ<sup>B</sup> p � f HGð Þþ cos θ<sup>B</sup> ð Þ� 1 � p f HGð Þ cos θ<sup>S</sup>

þ

� �<sup>2</sup> ( ): (11)

1 ffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>a</sup><sup>2</sup> <sup>þ</sup> <sup>1</sup> <sup>p</sup>

� �

1

This conversion formula was derived using probability theory [20].

where a is the bias coefficient that can be selected between 0 and 1. Once a biased angle θ<sup>B</sup> is randomly selected, away from the direction of the center of the

rotated around v^ by an angle ϕ randomly sampled from a uniform distribution in the range from 0 to 2π. These parameters are defined in the same manner as those used in the biased distribution presented in Section 3. The only difference is that the domain of cosð Þ θ<sup>B</sup> is restricted to a maximum deviation from the biased angle by 90°. This ensures that there would not be packets with very large likelihood ratio that could reduce the efficiency of our importance sampling. The likelihood ratio of the photon packet that uses the biased probability density function in Eq. (9) is

<sup>3</sup>=<sup>2</sup> , cosð Þ θ<sup>B</sup> ∈ ½ � 0, 1

, the provisional biased scattering direction u^<sup>0</sup> is

� � <sup>1</sup> <sup>þ</sup> <sup>a</sup><sup>2</sup> � <sup>2</sup><sup>a</sup> cosð Þ <sup>θ</sup><sup>B</sup>

1 þ g<sup>2</sup> � 2g cosð Þ θ<sup>S</sup> � �<sup>3</sup>=<sup>2</sup>

, (9)

,

(10)

: (12)

#### Figure 3.

Class I diffuse reflectance dependence on the distance from the center of the optical fiber for the simulation whose schematic is shown in Figure 4. The solid black line represents the result obtained with 2 <sup>10</sup><sup>5</sup> photon packets using the importance sampling technique presented in Section 3. The green long dashed line is the result obtained with 109 photon packets using standard MCML. The blue dots represent results obtained with 106 photon packets using MCML. The pink short dashed lines show the confidence interval of the importance sampling simulations with 2 105 photon packets that were estimated using a much larger ensemble of 64 105 simulations. Reprinted with permission from [14] © The Optical Society of America.

The second layer, extending from 2.32 to 2.42 mm from the tip of the fiber, has the same absorption and scattering coefficients as the first layer, but its refractive index is n = 1.33. The third layer, extending from 2.42 to 2.62 mm from the tip of the fiber, has the following parameters: μ<sup>a</sup> = 1.5 cm<sup>1</sup> , μ<sup>s</sup> = 30 cm<sup>1</sup> , and n = 1. After the third layer, the medium was assumed to be air: μ<sup>a</sup> = 0 cm<sup>1</sup> , μ<sup>s</sup> = 0 cm<sup>1</sup> , and n = 1. The anisotropy factor was assumed g = 0.9 for the three diffusive layers.

From Figure 3, we note an excellent correspondence between results obtained with our new importance sampling method and results obtained using MCML, that is, standard Monte Carlo simulations. However, our results were obtained in one-thousandth of the time required by the standard method.
