Two-step verification method for Monte Carlo codes in biomedical optics applications

Angelo Sassaroli; Federico Tommasi; Stefano Cavalieri; Lorenzo Fini; André Liemert; Alwin Kienle; Tiziano Binzoni; Fabrizio Martelli

doi:10.1117/1.JBO.27.8.083018

20 April 2022 Two-step verification method for Monte Carlo codes in biomedical optics applications

Angelo Sassaroli, Federico Tommasi, Stefano Cavalieri, Lorenzo Fini, André Liemert, Alwin Kienle, Tiziano Binzoni, Fabrizio Martelli

Author Affiliations +

Journal of Biomedical Optics, Vol. 27, Issue 8, 083018 (April 2022). https://doi.org/10.1117/1.JBO.27.8.083018

Abstract

Significance: Code verification is an unavoidable step prior to using a Monte Carlo (MC) code. Indeed, in biomedical optics, a widespread verification procedure for MC codes is still missing. Analytical benchmarks that can be easily used for the verification of different MC routines offer an important resource.

Aim: We aim to provide a two-step verification procedure for MC codes enabling the two main tasks of an MC simulator: (1) the generation of photons’ trajectories and (2) the intersections of trajectories with boundaries separating the regions with different optical properties. The proposed method is purely based on elementary analytical benchmarks, therefore, the correctness of an MC code can be assessed with a one-sample t-test.

Approach: The two-step verification is based on the following two analytical benchmarks: (1) the exact analytical formulas for the statistical moments of the spatial coordinates where the scattering events occur in an infinite medium and (2) the exact invariant solutions of the radiative transfer equation for radiance, fluence rate, and mean path length in media subjected to a Lambertian illumination.

Results: We carried out a wide set of comparisons between MC results and the two analytical benchmarks for a wide range of optical properties (from non-scattering to highly scattering media, with different types of scattering functions) in an infinite non-absorbing medium (step 1) and in a non-absorbing slab (step 2). The deviations between MC results and exact analytical values are usually within two standard errors (i.e., t-tests not rejected at a 5% level of significance). The comparisons show that the accuracy of the verification increases with the number of simulated trajectories so that, in principle, an arbitrary accuracy can be obtained.

Conclusions: Given the simplicity of the verification method proposed, we envision that it can be widely used in the field of biomedical optics.

1. Introduction

The verification of a Monte Carlo (MC) code is an important aspect of the whole process of confirmation that assures the scientific community of the reliability of its results.¹ In the process of confirmation, we can distinguish between a verification phase and a validation phase. The verification of an MC code is typically done by comparisons between its results and those obtained either with analytical benchmarks,¹^–³ or more often with previously verified MC codes. In contrast, the validation of an MC code is done by comparisons between its results and those obtained with experiments.⁴^,⁵ In this work, we are concerned only with the verification of an MC code. In the last decades, the modeling of heterogeneous tissue structures in MC codes for photon transport has required the development of algorithms with increasing complexity.⁶^–¹⁷ Therefore, the need for a thorough verification procedure has become more and more urgent.

As a matter of fact, most of the MC codes developed in biomedical optics have been verified partially or exclusively by means of comparisons with previously verified codes.⁶^–¹⁸ In the cross-verification procedures between MC results,⁹^–¹⁶^,¹⁸ the Monte Carlo modeling of light transport in multi-layered tissues (MCML) open-source code developed in the early nineties¹⁹ has been largely used as the standard reference. Due to this fact, the results of this code developed for a multi-layered medium can be considered as a numerical benchmark for photon migration through layered media. This special role played by MCML in biomedical applications is not observed for other MC codes.

One drawback of using verified MC codes to generate reference data is the limited accuracy that is achievable with a computer simulation. However, this drawback has not actually restricted the use of the MC method. In recent years, the facilitated access to open-source platforms of MC codes has made MC results easily available to a wider audience, and the practice of using verified MC codes has become the widest verification method used in biomedical applications.

It is also important to stress that in biomedical optics there is limited use of exact solutions of the radiative transfer equation (RTE) for the verification of MC codes. Indeed, only a few examples can be found for this kind of verification.¹²^,¹⁶^,¹⁹ This fact is related to the intrinsic complexity of both older and newly available solutions of the RTE,²⁰^–²³ which are not in closed-form and require some numerical evaluation. The verification of MCML, and other MC codes, has been largely based on the RTE solutions tabulated by van de Hulst,²⁴^,²⁵ related to a scattering slab, and on the RTE solutions for a semi-infinite medium tabulated by Giovanelli.²⁶ The extensive use of these historical results from Giovanelli and Van de Hulst, even though affected by a limited accuracy, provides another clear indicator of the difficulties to use other benchmarks based on more complex newly available solutions of the RTE.²⁰^–²³

The intent of this work is to propose a verification method purely based on the exact analytical solutions of the RTE. The selected benchmarks have the common characteristic to be extremely simple in terms of implementation so that the computational burden found with newly available RTE solutions²⁰^–²³ is avoided. Therefore, their use is also open to the non-expert users of complex computational methods since their implementation is straightforward. The method is framed in two steps which focus on two different aspects of photon migration in scattering media: (a) propagation through homogeneous domains, where the statistical rules valid in an infinite medium for the extraction of photons’ trajectories are used²^,²⁷^,²⁸ and (b) propagation through regions with different optical properties, where the effects of boundaries must be accounted for.¹⁹^,²⁷^,²⁸ Accordingly to this vision any verification method of an MC code can be in principle divided into (a) verification in an infinite medium where the basic routines/algorithms to extract the photons’ trajectories can be tested by using very intrinsic analytical benchmarks of photons transport and (b) verification in finite media, or in media with a least one finite dimension, where the effects of boundaries can be tested by means of specific benchmarks.

Thus, as the first step, we propose to verify an MC code by using statistical formulas for the first and the second moments of the spatial coordinates where the different scattering orders occur in an infinite non-absorbing medium.² For the second step of the verification method, we propose to use the invariant solutions for the radiance, $I$ , the fluence rate, $Φ$ , and mean total path length, $⟨ L ⟩$ , in bounded non-absorbing media subjected to Lambertian illumination.³^,²⁹^–³¹ With this two-step procedure, all the main quantities involved in an MC simulation can be verified with a high level of accuracy. Finally, it can be noted that this procedure is suitable to be used during the development of an MC code, i.e., from its very beginning till the finalization of the code. This aspect is not secondary in our proposed method and has an important didactic value: instead of verifying a code at the end of its implementation, a verification “in progress” has the advantage to check the different algorithms and routines during their construction with the advantage that the effects of multiple bugs can be better detected.

It is worth noting that the proposed method largely extends the verification methods previously published in Refs. 2 and 3. The benchmarks used in Ref. 2 for isotropic scattering can now be applied for any scattering order. Compared to Ref. 3, the benchmarks used include fundamental radiometric quantities such as radiance and fluence rate. This latter extension is of fundamental interest because the radiance is at the origin of all the radiometric quantities used in radiative transfer, therefore a verification based on radiance is quite relevant. Moreover, the fluence rate represents a fundamental quantity in the latest generation of MC codes describing propagation in complex geometries,⁶^–¹⁴^,¹⁶^,¹⁷ therefore, a verification method based on fluence rate is of the utmost importance.

In Sec. 2, the analytical benchmarks used in the proposed verification method are described. In the same section, the proposed verification method is also described. In Sec. 3, the results obtained with the two-step verification method are presented. In Sec. 4, the results described in Sec. 3 are briefly discussed together with the future perspectives for the application of this method.

2. Theory and Methods

This section describes the two kinds of benchmarks used in the verification procedure proposed in this work together with a brief description of the method.

2.1.

Analytical Benchmarks for Light Propagation through an Infinite Medium: Statistical Moments of the Coordinates of the Scattering Events

In this section, the relationships existing between the statistical moments of the coordinates where scattering events occur and the optical properties are presented for the case of an infinite non-absorbing medium. It is at first considered the general case of rotationally symmetric scattering phase functions, i.e., scattering functions which can be represented as $p (\hat{s} \cdot {\hat{s}}^{'})$ , with $\hat{s}$ direction of the incident radiation and ${\hat{s}}^{'}$ direction of the scattered radiation (Sec. 2.1.1). Then, it is considered the special case of isotropic scattering (Appendix). The benchmarks presented in this section are a set of analytical solutions of immediate implementation. The proof of their validity can be found in the related cited literature (see below) and in the Appendix included in this work.

2.1.1.

Rotationally symmetric scattering phase functions

Let’s consider an infinite homogeneous non-absorbing medium, where at the origin of a Cartesian reference system we have a pencil light source in the $z$ direction. The optical properties of the medium are described by the scattering coefficient, $μ_{s}$ , and the scattering function, $p (θ)$ , with $θ$ the scattering angle.²⁷^,²⁸ Although we consider a non-absorbing medium, it is important to note that the inclusion of absorption is straightforward by means of the microscopic Beer–Lambert law (mBLL).²⁸ In what follows the average values of the first and second moments of the coordinates where the different scattering orders occur (up to the fourth-order) are reported.² The different scattering orders will be denoted with a subscript giving the number of the scattering order. Thus, the mean values of $x_{1}$ , $y_{1}$ , and $z_{1}$ (first moments) where the first scattering event occurs are given by²

Eq. (1)

⟨ x_{1} ⟩ = ⟨ y_{1} ⟩ = 0,

Eq. (2)

⟨ z_{1} ⟩ = \frac{1}{μ_{s}},

and the relative mean square values

⟨ x_{1}^{2} ⟩

,

⟨ y_{1}^{2} ⟩

, and

⟨ z_{1}^{2} ⟩

(second moments) are given by²

Eq. (3)

⟨ x_{1}^{2} ⟩ = ⟨ y_{1}^{2} ⟩ = 0,

Eq. (4)

⟨ z_{1}^{2} ⟩ = \frac{2}{μ_{s}^{2}} .

From the above relations, the values of the mean square distance from the

z

axis,

⟨ ρ_{1}^{2} ⟩

, and the mean square distance from the source,

⟨ d_{1}^{2} ⟩

, are obtained as²

Eq. (5)

⟨ ρ_{1}^{2} ⟩ = ⟨ x_{1}^{2} + y_{1}^{2} ⟩ = ⟨ x_{1}^{2} ⟩ + ⟨ y_{1}^{2} ⟩ = 0,

Eq. (6)

⟨ d_{1}^{2} ⟩ = ⟨ x_{1}^{2} + y_{1}^{2} + z_{1}^{2} ⟩ = ⟨ x_{1}^{2} ⟩ + ⟨ y_{1}^{2} ⟩ + ⟨ z_{1}^{2} ⟩ = \frac{2}{μ_{s}^{2}} .

For the second scattering order, we have²

Eq. (7)

⟨ x_{2} ⟩ = ⟨ y_{2} ⟩ = 0,

Eq. (8)

⟨ z_{2} ⟩ = \frac{1 + g}{μ_{s}},

Eq. (9)

⟨ x_{2}^{2} ⟩ = ⟨ y_{2}^{2} ⟩ = \frac{1 - g_{2}}{μ_{s}^{2}},

Eq. (10)

⟨ z_{2}^{2} ⟩ = \frac{2 (1 + g + g_{2})}{μ_{s}^{2}},

where

g = ⟨ \cos θ ⟩

and

g_{2} = ⟨ \cos^{2} θ ⟩

are the first and second moments of the cosine of the scattering angle, respectively. The parameter

g

is usually denoted asymmetry factor of the phase scattering function.

From Eqs. (9) and (10) we have²

Eq. (11)

⟨ ρ_{2}^{2} ⟩ = ⟨ x_{2}^{2} + y_{2}^{2} ⟩ = ⟨ x_{2}^{2} ⟩ + ⟨ y_{2}^{2} ⟩ = \frac{2 (1 - g_{2})}{μ_{s}^{2}},

Eq. (12)

⟨ d_{2}^{2} ⟩ = ⟨ x_{2}^{2} + y_{2}^{2} + z_{2}^{2} ⟩ = ⟨ x_{2}^{2} ⟩ + ⟨ y_{2}^{2} ⟩ + ⟨ z_{2}^{2} ⟩ = \frac{2 (2 + g)}{μ_{s}^{2}} .

The mean coordinates of the point at which the third-order of scattering occurs are²

Eq. (13)

⟨ x_{3} ⟩ = ⟨ y_{3} ⟩ = 0,

Eq. (14)

⟨ z_{3} ⟩ = \frac{1 + g + g^{2}}{μ_{s}},

and also

Eq. (15)

⟨ x_{3}^{2} ⟩ = ⟨ y_{3}^{2} ⟩ = \frac{\frac{3}{2} + g - g g_{2} - \frac{3}{2} g_{2}^{2}}{μ_{s}^{2}},

Eq. (16)

⟨ z_{3}^{2} ⟩ = \frac{3 + 2 g + 2 g^{2} + 2 g g_{2} + 3 g_{2}^{2}}{μ_{s}^{2}} .

Similarly to the previous scattering orders, we also have

Eq. (17)

⟨ ρ_{3}^{2} ⟩ = ⟨ x_{3}^{2} + y_{3}^{2} ⟩ = ⟨ x_{3}^{2} ⟩ + ⟨ y_{3}^{2} ⟩ = \frac{3 + 2 g - 2 g g_{2} - 3 g_{2}^{2}}{μ_{s}^{2}},

Eq. (18)

⟨ d_{3}^{2} ⟩ = ⟨ x_{3}^{2} + y_{3}^{2} + z_{3}^{2} ⟩ = ⟨ x_{3}^{2} ⟩ + ⟨ y_{3}^{2} ⟩ + ⟨ z_{3}^{2} ⟩ = \frac{2 (3 + 2 g + g^{2})}{μ_{s}^{2}} .

For the fourth-order, we have²

Eq. (19)

⟨ x_{4} ⟩ = ⟨ y_{4} ⟩ = 0,

Eq. (20)

⟨ z_{4} ⟩ = \frac{1 + g + g^{2} + g^{3}}{μ_{s}},

Eq. (21)

⟨ x_{4}^{2} ⟩ = ⟨ y_{4}^{2} ⟩ = \frac{\frac{9}{4} + \frac{3}{2} g - \frac{3}{2} g g_{2}^{2} + g^{2} - g^{2} g_{2} - \frac{3}{4} g_{2} + \frac{3}{4} g_{2}^{2} - \frac{9}{4} g_{2}^{3}}{μ_{s}^{2}},

Eq. (22)

⟨ z_{4}^{2} ⟩ = \frac{\frac{7}{2} + 3 g (1 + g_{2}^{2}) + 2 g^{2} (1 + g_{2}) + 2 g^{3} + \frac{3}{2} g_{2} - \frac{3}{2} g_{2}^{2} + \frac{9}{2} g_{2}^{3}}{μ_{s}^{2}} .

Eq. (23)

⟨ ρ_{4}^{2} ⟩ = ⟨ x_{4}^{2} + y_{4}^{2} ⟩ = ⟨ x_{4}^{2} ⟩ + ⟨ y_{4}^{2} ⟩ = \frac{\frac{9}{2} + 3 g - 3 g g_{2}^{2} + 2 g^{2} - 2 g^{2} g_{2} + - \frac{3}{2} g_{2} + \frac{3}{2} g_{2}^{2} - \frac{9}{2} g_{2}^{3}}{μ_{s}^{2}},

Eq. (24)

⟨ d_{4}^{2} ⟩ = ⟨ x_{4}^{2} + y_{4}^{2} + z_{4}^{2} ⟩ = ⟨ x_{4}^{2} ⟩ + ⟨ y_{4}^{2} ⟩ + ⟨ z_{4}^{2} ⟩ = \frac{2 (4 + 3 g + 2 g^{2} + g^{3})}{μ_{s}^{2}} .

The calculations can be in principle iterated for the higher scattering orders. Given the cylindrical symmetry, with respect to the $z$ axis, we have that $⟨ x_{k} ⟩$ and $⟨ y_{k} ⟩$ are zero for any $k$ . While for $⟨ z_{k} ⟩$ , the following relation is valid, which is given by²

Eq. (25)

⟨ z_{k} ⟩ = \frac{1 - g^{k}}{μ_{s} (1 - g)} .

It is worth noting that the above equation for

k \to \infty

returns the classical transport mean free path,

\frac{1}{μ_{s} (1 - g)}

.²⁸ For the higher scattering orders it is also possible to extrapolate the following recurring relation for

⟨ d_{k}^{2} ⟩

, which is given by²

Eq. (26)

⟨ d_{k}^{2} ⟩ = ⟨ x_{k}^{2} + y_{k}^{2} + z_{k}^{2} ⟩ = 2 \frac{k - (k + 1) g + g^{k + 1}}{μ_{s}^{2} {(1 - g)}^{2}} .

A full proof of Eq. (26), derived from the RTE, was obtained by Liemert et al.³² for arbitrary rotationally symmetric scattering phase functions (e.g., the Henyey–Greenstein (HG) phase function).²⁸ In Ref. 32, this expression is also generalized for an absorbing medium. The importance of this benchmark is that it can be used for any scattering order and is usually employed in biomedical optics for most of the scattering functions.

Statistical relationships can also be obtained for the optical path lengths of the propagated photons. The statistical moment of order $m$ of the path length traveled at the $k$ ’th scattering event is given by²

Eq. (27)

⟨ l_{k}^{m} ⟩ = \frac{k (k + 1) \dots (k + m - 1)}{μ_{s}^{m}} .

2.1.2.

Isotropic scattering phase functions

For the case of isotropic scattering, it is possible to extrapolate the following recurring relations for the second moment for any scattering order $k > 0$ as given by

Eq. (28)

⟨ x_{k}^{2} ⟩ = ⟨ y_{k}^{2} ⟩ = \frac{2}{3 μ_{s}^{2}} (k - 1),

Eq. (29)

⟨ z_{k}^{2} ⟩ = \frac{2}{3 μ_{s}^{2}} (k + 2),

Eq. (30)

⟨ ρ_{k}^{2} ⟩ = ⟨ x_{k}^{2} + y_{k}^{2} ⟩ = ⟨ x_{k}^{2} ⟩ + ⟨ y_{k}^{2} ⟩ = \frac{4}{3 μ_{s}^{2}} (k - 1),

Eq. (31)

⟨ d_{k}^{2} ⟩ = ⟨ x_{k}^{2} + y_{k}^{2} + z_{k}^{2} ⟩ = ⟨ x_{k}^{2} ⟩ + ⟨ y_{k}^{2} ⟩ + ⟨ z_{k}^{2} ⟩ = \frac{2 k}{μ_{s}^{2}} .

The proof of the validity of the above equations is treated in the Appendix where exact analytical equations for the moments of the coordinates where scattering events occur are obtained. We notice that, compared to the previous benchmarks (Sec. 2.1.1), the calculated moments for isotropic scattering do not depend on

g

and

g_{2}

.

The presented benchmarks provide an overview of the “cloud” of photons migrating in an infinite medium. The first moments yield the coordinates of its barycenter, while the second central moments represent its width along the three axes. They are strictly related to the scattering function and the statistical law for scattering interactions. The features of the scattering function that affect the statistical moments are: $g = ⟨ \cos θ ⟩$ and $g_{2} = ⟨ \cos^{2} θ ⟩$ . Thus, these benchmarks are suitable to verify if an MC code correctly extracts the photons’ trajectories and consequently the scattering function is correctly implemented inside the code.

2.2.

Analytical Benchmarks for Boundary Effects: Invariant Solutions for Radiance, Fluence Rate, and Total Mean Path Length

Photon migration in a finite medium is also characterized by boundary effects, i.e., reflection, refraction, and intersection. Besides the boundary with the outer medium, internal boundaries are also used to enclose regions with different optical properties. The correct evaluation of boundary effects is another important task of an MC code. To verify the reliability of an MC code to simulate boundary effects, we propose to use a set of invariant exact solutions of the RTE for the radiance, fluence rate, and for the partial and total path lengths that are obtained when a non-absorbing medium is subjected to a Lambertian illumination.²⁹^,³⁰^,³³^,³⁴ These solutions depend on the distribution of the refractive index inside the medium relative to the refractive index of the outer medium.

Let’s consider a non-absorbing inhomogeneous medium of volume $V$ (with at least one finite dimension) without internal sources, delimited by a smooth convex surface $Σ$ . The medium is composed of a number $N$ of discrete sub-volumes $V_{j}$ of refractive index $n_{j}$ and with no restrictions on the scattering properties inside each $V_{j}$ . The surfaces enclosing each sub-volume are assumed to be smooth so that Snell’s and Fresnel’s law can be applied. The refractive index of the external medium is denoted with $n_{e}$ and the surface $Σ$ is illuminated by a continuous wave Lambertian radiation of intensity $I_{0} (W m^{- 2} {sr}^{- 1})$ , i.e., the source term is thus can be expressed in terms of the distribution of radiance on the external surface $Σ$ , which is given by

Eq. (32)

I_{eSource} ({\vec{r}}_{Σ}, {\hat{s}}_{e}) = I_{0} [W m^{- 2} {sr}^{- 1}] \forall {\vec{r}}_{Σ} \in Σ and \forall {\hat{s}}_{e} inwardly directed to Σ .

In this case, the solution for the radiance inside

V_{j}

is²⁹

Eq. (33)

I_{j} (\vec{r}, \hat{s}) = {(\frac{n_{j}}{n_{e}})}^{2} I_{0} [W m^{- 2} {sr}^{- 1}] \forall \hat{s} and \forall \vec{r} \in V_{j},

where

I_{j} (\vec{r}, \hat{s})

is a function of the refractive indices

n_{j}

, internal to

V_{j}

, and

n_{e}

, external to

V

; however, we notice that this solution does not depend on the refractive index of the remaining sub-volumes. Also, in Eq. (33)

\hat{s}

denotes the direction vector, while

\vec{r}

denotes the position vector internal to the medium. Consequently, the solution for the fluence inside

V_{j}

is expressed as²⁹

Eq. (34)

Φ_{j} (\vec{r}) = \int_{4 π} I_{j} (\vec{r}, \hat{s}) d \hat{s} = 4 π {(\frac{n_{j}}{n_{e}})}^{2} I_{0}, [W m^{- 2}], \forall \vec{r} \in V_{j} .

Finally, the average internal path length in each sub-volume

V_{j}

results in²⁹

Eq. (35)

⟨ L_{j} ⟩ = 4 {(\frac{n_{j}}{n_{e}})}^{2} \frac{V_{j}}{Σ} .

The average total path length

⟨ L ⟩

inside the total volume

V

can also be evaluated by simply summing up all the contributions of the average internal path lengths

⟨ L_{j} ⟩

of Eq. (35), i.e.,

Eq. (36)

⟨ L ⟩ = \sum_{j = 1}^{N} ⟨ L_{j} ⟩ = 4 \sum_{j = 1}^{N} {(\frac{n_{j}}{n_{e}})}^{2} \frac{V_{j}}{Σ} .

From Eqs. (33) and (34), we clearly see that the solutions for the internal radiance and the fluence rate are invariant both with respect to the geometry of the medium (provided that the external surface is convex and smooth) and with respect to the distribution of the scattering properties (scattering coefficient, scattering function, and homogeneity). Similarly, the expressions for the average internal path lengths

⟨ L_{j} ⟩

and the total path length

⟨ L ⟩

are invariant with respect to the scattering properties of the medium, however, they depend (as it is expected) on the volumes of the sub-regions and the whole external surface

Σ

. One exception is the case of

μ_{s j} (\vec{r}) = 0

whenever the geometry (e.g., sphere or slab) of the medium allows for a regime of trapped trajectories.²⁹ These properties imply that these benchmarks can be used for any scattering coefficient of the medium except, with

μ_{s j} (\vec{r}) = 0

, for all those geometries where trapped trajectories can be established.

In general, for a non-scattering medium [ $μ_{s j} (\vec{r}) = 0$ ] with $n_{j} \leq n_{e}$ , the above solutions are still valid. While, with $μ_{s j} (\vec{r}) = 0$ , for geometries like a sphere or slab with $n_{j} > n_{e}$ ,²⁹^,³⁵ where photons can be trapped inside the volume, the same RTE solutions found for scattering media cannot be used. For example, the RTE solution for a layered non-scattering slab (which will be used in the results section) should account for the regime of trapped photons. Accordingly, the RTE solution for the internal radiance inside a layered non-scattering slab is²⁹

Eq. (37)

I_{j} (\vec{r}, \hat{s}) = {(\frac{n_{j}}{n_{e}})}^{2} I_{0} \forall \vec{r} \in V_{j} \forall \hat{s} | | \hat{s} \cdot \hat{q} | \geq \cos θ_{j Ma x} I_{j} (\vec{r}, \hat{s}) = 0 \forall \vec{r} \in V_{j} \forall \hat{s} | | \hat{s} \cdot \hat{q} | < \cos θ_{j Max},

where

I_{0}

is the radiance on the external surface

Σ

,

\hat{q}

is the unit vector perpendicular to the slab inwardly directed,

θ_{j Max}

is the maximum entrance angle in the medium for having radiation in the

j

’th layer. Consequently, the solution for the fluence rate in a layered slab is

Eq. (38)

Φ_{j} (\vec{r}) = 4 π {(\frac{n_{j}}{n_{e}})}^{2} I_{0} [1 - \cos (θ_{j Max})], [W m^{- 2}] \forall \vec{r} \in V_{j} .

Finally, the average internal path length

⟨ L_{j} ⟩

spent in the

j

’th layer of the slab is

Eq. (39)

⟨ L_{j} ⟩ = 2 s_{j} {(\frac{n_{j}}{n_{e}})}^{2} [1 - \cos (θ_{j Max})],

where

s_{j}

is the thickness of the

j

’th layer of the slab. For some guidelines to calculate the angle

θ_{j Max}

, we refer to the work of Martelli et al.²⁹ It may be worth noting that for

n_{j} > n_{e}

in the above solutions we have a discontinuity of the radiance between the case

μ_{s} = 0

[Eq. (37)] and the case

μ_{s} \neq 0

[Eq. (33)].²⁹ When

n_{j} > n_{e}

, a discontinuity is also observed for the fluence rate,

Φ_{j}

and the partial path length

⟨ L_{j} ⟩

as can be similarly noted by the above equations. ²⁹^,³⁵ Moreover, when

n_{j} > n_{e}

and

μ_{s} = 0

, from Eq. (37) the discontinuity of the radiance can be observed for the angle

θ_{j Max}

. At such value of the angle the radiance

I_{j}

switches from the value

(\frac{n_{j}}{n_{e}})^{2} I_{0}

to zero.

The presented benchmarks provide an overview of invariant solutions of the RTE in non-absorbing media that are sensitive to the effects of boundary conditions between the different regions of the medium and also between the medium and the external region. It must be noted that the independence of the presented solutions from the scattering properties of the medium implies that they cannot be used to test the correctness of the phase function implemented in the code. However, this is done with the previous benchmark of Sec. 2.1. Combined together, the benchmarks of Secs. 2.1 and 2.2 cover the main characteristics of photon migration according to the RTE.

2.3.

Proposed Method

The verification method here proposed is divided into two steps and exploits the following strategy: (1) all the routines of the code involved in the extraction of the trajectories, based also on the scattering function are verified by a direct comparison of MC results in an infinite medium with the benchmarks as in Sec. 2.1 and (2) the routines of the MC code related to the interactions of photons’ trajectories with boundaries are verified by a direct comparison of the MC results in a slab geometry with the invariant solutions for radiance, fluence rate, and mean path length as in Sec. 2.2.

Concisely, the first step is devised to test the phase function implemented in the MC code, the second one to test the intersection of photon’s trajectories with boundaries, including the correct implementation of the reflection and refraction laws. The second step includes the boundary with the external medium and those between regions of the medium with different optical properties. We notice that the correct calculation of the intersections with boundaries is at the core of partial path length estimation, which is one fundamental task of an MC code. In fact, the absence of boundaries in the first step is ideal for testing the phase function. In contrast, the invariant solutions of RTE used in the second step (which do not depend on the features of the phase function) are ideal for testing the intersection with boundaries.

We believe that the two-step verification increases the sensitivity to detect errors in an MC code.

3. Results

In this section, the results obtained using the proposed method in the two steps of the verification are presented. The results obtained in the first step of the verification are described in Sec. 3.1. While the results obtained in the second step of the verification are described in Sec. 3.2.

3.1.

First step of the Verification (Sec. 2.1)

In this section, the moments of the coordinates of the scattering points calculated with MC simulations have been compared with the analytical benchmarks of Sec. 2.1. All the comparisons are carried out in an infinite non-absorbing medium with a unitary scattering coefficient $μ_{s} = 1 {mm}^{- 1}$ where a pencil beam source is injected at the origin of a Cartesian reference system along the $z$ direction. Three scattering functions have been selected: the HG model with $g = 0$ and 0.9,²⁸ and a Rayleigh scattering function ( $g = 0$ ).²⁸ This choice is motivated on this ground: for these scattering functions the moments $g$ and $g_{2}$ are known exactly without resorting to numerical evaluations. This is not true for phase scattering functions derived from Mie theory.²⁸

The MC code here subjected to the verification procedure is a program used for simulating light propagation in tissue optics. The core of the program, developed during the period 1980 to 1999,²⁸^,³⁶^–³⁹ generates a large number of photons’ trajectories for different scattering coefficients and refractive indices. In Tables 1 Table 2–3 the relative errors (ratio between standard error and average value) of the MC results obtained for the moments of Sec. 2.1 for the first four scattering orders ( $k \in {1, \dots, 4}$ ) are shown for three values of the number of simulated trajectories $N \in {10^{6}, 10^{8}, 10^{10}}$ . We notice that the relative error decreases by one order of magnitude as the value of $N$ increases by two orders of magnitude. This behavior expresses exactly the expected distributions for the calculated quantities and highlights the fact that the precision of the calculation is only limited by the numerical accuracy of the processor and by the number of simulated trajectories, as is expected from the statistical theory. In these tables, the data for $⟨ x_{1} ⟩$ , $⟨ y_{1} ⟩$ , $⟨ x_{1}^{2} ⟩$ , $⟨ y_{1}^{2} ⟩$ , and $⟨ ρ_{1}^{2} ⟩$ have not been reported since their values are null and the relative uncertainty cannot be calculated.

Table 1

Relative uncertainty, σ⟨…⟩MC⟨…⟩MC, of the statistical moments of the coordinates of scattering events calculated with MC simulations for a non-absorbing infinite medium with scattering coefficient μs=1 mm−1, isotropic scattering function obtained with the HG model (g=0, g2=1/3), and three values N of generated trajectories.

σ⟨…⟩MC⟨…⟩MC	N=106				N=108				N=1010
k	1	2	3	4	1	2	3	4	1	2	3	4
$⟨ z_{k} ⟩$	$1.000 \cdot 10^{- 3}$	$1.291 \cdot 10^{- 3}$	$1.529 \cdot 10^{- 3}$	$1.732 \cdot 10^{- 3}$	$1.000 \cdot 10^{- 4}$	$1.291 \cdot 10^{- 4}$	$1.528 \cdot 10^{- 4}$	$1.732 \cdot 10^{- 4}$	$1.000 \cdot 10^{- 5}$	$1.291 \cdot 10^{- 5}$	$1.528 \cdot 10^{- 5}$	$1.732 \cdot 10^{- 5}$
$⟨ x_{k}^{2} ⟩$		$3.125 \cdot 10^{- 3}$	$2.436 \cdot 10^{- 3}$	$2.152 \cdot 10^{- 3}$		$3.133 \cdot 10^{- 4}$	$2.429 \cdot 10^{- 4}$	$2.145 \cdot 10^{- 4}$		$3.130 \cdot 10^{- 5}$	$2.429 \cdot 10^{- 5}$	$2.145 \cdot 10^{- 5}$
$⟨ y_{k}^{2} ⟩$		$3.134 \cdot 10^{- 3}$	$2.428 \cdot 10^{- 3}$	$2.141 \cdot 10^{- 3}$		$3.132 \cdot 10^{- 4}$	$2.429 \cdot 10^{- 4}$	$2.145 \cdot 10^{- 4}$		$3.131 \cdot 10^{- 5}$	$2.429 \cdot 10^{- 5}$	$2.145 \cdot 10^{- 5}$
$⟨ z_{k}^{2} ⟩$	$2.229 \cdot 10^{- 3}$	$2.042 \cdot 10^{- 3}$	$1.927 \cdot 10^{- 3}$	$1.845 \cdot 10^{- 3}$	$2.236 \cdot 10^{- 4}$	$2.043 \cdot 10^{- 4}$	$1.925 \cdot 10^{- 4}$	$1.845 \cdot 10^{- 4}$	$2.236 \cdot 10^{- 5}$	$2.043 \cdot 10^{- 5}$	$1.925 \cdot 10^{- 5}$	$1.844 \cdot 10^{- 5}$
$⟨ ρ_{k}^{2} ⟩$		$2.486 \cdot 10^{- 3}$	$1.899 \cdot 10^{- 3}$	$1.655 \cdot 10^{- 3}$		$2.492 \cdot 10^{- 4}$	$1.898 \cdot 10^{- 4}$	$1.653 \cdot 10^{- 4}$		$2.490 \cdot 10^{- 5}$	$1.897 \cdot 10^{- 5}$	$1.653 \cdot 10^{- 5}$
$⟨ d_{k}^{2} ⟩$	$2.229 \cdot 10^{- 3}$	$1.682 \cdot 10^{- 3}$	$1.455 \cdot 10^{- 3}$	$1.323 \cdot 10^{- 3}$	$2.236 \cdot 10^{- 4}$	$. 684 \cdot 10^{- 4}$	$1.453 \cdot 10^{- 4}$	$1.323 \cdot 10^{- 4}$	$2.236 \cdot 10^{- 5}$	$1.683 \cdot 10^{- 5}$	$1.453 \cdot 10^{- 5}$	$1.323 \cdot 10^{- 5}$
$⟨ l_{k} ⟩$	$1.000 \cdot 10^{- 3}$	$7.077 \cdot 10^{- 4}$	$5.776 \cdot 10^{- 4}$	$4.999 \cdot 10^{- 4}$	$1.000 \cdot 10^{- 4}$	$7.071 \cdot 10^{- 5}$	$5.773 \cdot 10^{- 5}$	$5.001 \cdot 10^{- 5}$	$1.000 \cdot 10^{- 5}$	$7.071 \cdot 10^{- 6}$	$5.774 \cdot 10^{- 6}$	$5.000 \cdot 10^{- 6}$

Table 2

Relative uncertainty, σ⟨…⟩MC⟨…⟩MC, of the statistical moments of the coordinates of scattering events calculated with MC simulations for a non-absorbing infinite medium with scattering coefficient μs=1 mm−1, an anisotropic scattering function obtained with the HG model with g=0.9 (g2=2.62/3), and three values N of generated trajectories.

σ⟨…⟩MC⟨…⟩MC	N=106				N=108				N=1010
k	1	2	3	4	1	2	3	4	1	2	3	4
$⟨ z_{k} ⟩$	$1.000 \cdot 10^{- 3}$	$7.330 \cdot 10^{- 4}$	$6.331 \cdot 10^{- 4}$	$5.884 \cdot 10^{- 4}$	$1 \cdot 10^{- 4}$	$7.324 \cdot 10^{- 6}$	$6.322 \cdot 10^{- 6}$	$5.880 \cdot 10^{- 5}$	$1.000 \cdot 10^{- 5}$	$7.324 \cdot 10^{- 6}$	$6.323 \cdot 10^{- 6}$	$5.880 \cdot 10^{- 6}$
$⟨ x_{k}^{2} ⟩$		$6.019 \cdot 10^{- 3}$	$3.795 \cdot 10^{- 3}$	$2.957 \cdot 10^{- 3}$		$6.009 \cdot 10^{- 4}$	$3.847 \cdot 10^{- 4}$	$2.986 \cdot 10^{- 4}$		$6.009 \cdot 10^{- 5}$	$3.846 \cdot 10^{- 5}$	$2.986 \cdot 10^{- 5}$
$⟨ y_{k}^{2} ⟩$		$5.931 \cdot 10^{- 3}$	$3.801 \cdot 10^{- 3}$	$2.967 \cdot 10^{- 3}$		$5.993 \cdot 10^{- 4}$	$3.845 \cdot 10^{- 4}$	$2.983 \cdot 10^{- 4}$		$6.007 \cdot 10^{- 5}$	$3.845 \cdot 10^{- 5}$	$2.986 \cdot 10^{- 5}$
$⟨ z_{k}^{2} ⟩$	$2.229 \cdot 10^{- 3}$	$1.564 \cdot 10^{- 3}$	$1.294 \cdot 10^{- 3}$	$1.148 \cdot 10^{- 3}$	$2.236 \cdot 10^{- 4}$	$1.564 \cdot 10^{- 4}$	$1.293 \cdot 10^{- 4}$	$1.148 \cdot 10^{- 4}$	$2.236 \cdot 10^{- 5}$	$1.564 \cdot 10^{- 5}$	$1.293 \cdot 10^{- 5}$	$1.148 \cdot 10^{- 5}$
$⟨ ρ_{k}^{2} ⟩$		$4.838 \cdot 10^{- 3}$	$3.054 \cdot 10^{- 3}$	$2.352 \cdot 10^{- 3}$		$4.865 \cdot 10^{- 4}$	$3.086 \cdot 10^{- 4}$	$2.367 \cdot 10^{- 4}$		$4.871 \cdot 10^{- 5}$	$3.086 \cdot 10^{- 5}$	$2.369 \cdot 10^{- 5}$
$⟨ d_{k}^{2} ⟩$	$2.229 \cdot 10^{- 3}$	$1.528 \cdot 10^{- 3}$	$1.225 \cdot 10^{- 3}$	$1.048 \cdot 10^{- 3}$	$2.236 \cdot 10^{- 4}$	$1.538 \cdot 10^{- 4}$	$1.244 \cdot 10^{- 4}$	$1.076 \cdot 10^{- 4}$	$2.236 \cdot 10^{- 5}$	$1.538 \cdot 10^{- 5}$	$1.244 \cdot 10^{- 5}$	$1.076 \cdot 10^{- 5}$
$⟨ l_{k} ⟩$	$1.000 \cdot 10^{- 3}$	$7.077 \cdot 10^{- 3}$	$5.776 \cdot 10^{- 4}$	$4.999 \cdot 10^{- 4}$	$1.000 \cdot 10^{- 4}$	$7.071 \cdot 10^{- 5}$	$5.773 \cdot 10^{- 5}$	$5.001 \cdot 10^{- 5}$	$1.000 \cdot 10^{- 5}$	$7.071 \cdot 10^{- 6}$	$5.774 \cdot 10^{- 6}$	$5.000 \cdot 10^{- 6}$
$⟨ l_{k}^{2} ⟩$	$2.229 \cdot 10^{- 3}$	$1.528 \cdot 10^{- 3}$	$1.225 \cdot 10^{- 3}$	$1.048 \cdot 10^{- 3}$	$2.236 \cdot 10^{- 4}$	$1.528 \cdot 10^{- 4}$	$1.225 \cdot 10^{- 4}$	$1.049 \cdot 10^{- 4}$	$2.236 \cdot 10^{- 5}$	$1.528 \cdot 10^{- 5}$	$1.225 \cdot 10^{- 5}$	$1.049 \cdot 10^{- 5}$

Table 3

Relative uncertainty, σ⟨…⟩MC⟨…⟩MC, of the statistical moments of the coordinates of scattering events calculated with MC simulations for a non-absorbing infinite medium with scattering coefficient μs=1 mm−1, a Rayleigh scattering function (g=0, g2=0.4), and three values N of generated trajectories.

σ⟨…⟩MC⟨…⟩MC	N=106				N=108				N=1010
k	1	2	3	4	1	2	3	4	1	2	3	4
$⟨ z_{k} ⟩$	$1.000 \cdot 10^{- 3}$	$1.339 \cdot 10^{- 3}$	$1.573 \cdot 10^{- 3}$	$1.770 \cdot 10^{- 3}$	$1.000 \cdot 10^{- 4}$	$1.342 \cdot 10^{- 4}$	$1.575 \cdot 10^{- 4}$	$1.774 \cdot 10^{- 4}$	$1.000 \cdot 10^{- 5}$	$1.342 \cdot 10^{- 5}$	$1.575 \cdot 10^{- 5}$	$1.774 \cdot 10^{- 5}$
$⟨ x_{k}^{2} ⟩$		$3.222 \cdot 10^{- 3}$	$2.475 \cdot 10^{- 3}$	$2.185 \cdot 10^{- 3}$		$3.232 \cdot 10^{- 4}$	$2.490 \cdot 10^{- 4}$	$2.190 \cdot 10^{- 4}$		$3.229 \cdot 10^{- 5}$	$2.488 \cdot 10^{- 5}$	$2.189 \cdot 10^{- 5}$
$⟨ y_{k}^{2} ⟩$		$3.233 \cdot 10^{- 3}$	$2.493 \cdot 10^{- 3}$	$2.177 \cdot 10^{- 3}$		$3.231 \cdot 10^{- 4}$	$2.488 \cdot 10^{- 4}$	$2.189 \cdot 10^{- 4}$		$3.229 \cdot 10^{- 5}$	$2.488 \cdot 10^{- 5}$	$2.188 \cdot 10^{- 5}$
$⟨ z_{k}^{2} ⟩$	$2.229 \cdot 10^{- 3}$	$2.016 \cdot 10^{- 3}$	$1.912 \cdot 10^{- 3}$	$1.837 \cdot 10^{- 3}$	$2.236 \cdot 10^{- 4}$	$2.018 \cdot 10^{- 4}$	$1.912 \cdot 10^{- 4}$	$1.839 \cdot 10^{- 4}$	$2.236 \cdot 10^{- 5}$	$2.018 \cdot 10^{- 5}$	$1.912 \cdot 10^{- 5}$	$1.839 \cdot 10^{- 5}$
$⟨ ρ_{k}^{2} ⟩$		$2.567 \cdot 10^{- 3}$	$1.942 \cdot 10^{- 3}$	$1.686 \cdot 10^{- 3}$		$2.575 \cdot 10^{- 4}$	$1.948 \cdot 10^{- 4}$	$1.691 \cdot 10^{- 4}$		$2.573 \cdot 10^{- 5}$	$1.948 \cdot 10^{- 5}$	$1.691 \cdot 10^{- 5}$
$⟨ d_{k}^{2} ⟩$	$2.229 \cdot 10^{- 3}$	$1.702 \cdot 10^{- 3}$	$1.472 \cdot 10^{- 3}$	$1.339 \cdot 10^{- 3}$	$2.236 \cdot 10^{- 4}$	$1.703 \cdot 10^{- 4}$	$1.474 \cdot 10^{- 4}$	$1.343 \cdot 10^{- 4}$	$2.236 \cdot 10^{- 5}$	$1.703 \cdot 10^{- 5}$	$1.474 \cdot 10^{- 5}$	$1.343 \cdot 10^{- 5}$
$⟨ l_{k} ⟩$	$1.000 \cdot 10^{- 3}$	$7.077 \cdot 10^{- 3}$	$5.776 \cdot 10^{- 4}$	$4.999 \cdot 10^{- 4}$	$1.000 \cdot 10^{- 4}$	$7.071 \cdot 10^{- 5}$	$5.773 \cdot 10^{- 5}$	$5.001 \cdot 10^{- 5}$	$1.000 \cdot 10^{- 5}$	$7.071 \cdot 10^{- 6}$	$5.774 \cdot 10^{- 6}$	$5.000 \cdot 10^{- 6}$
$⟨ l_{k}^{2} ⟩$	$2.229 \cdot 10^{- 3}$	$1.528 \cdot 10^{- 3}$	$1.225 \cdot 10^{- 3}$	$1.048 \cdot 10^{- 3}$	$2.236 \cdot 10^{- 4}$	$1.528 \cdot 10^{- 4}$	$1.225 \cdot 10^{- 4}$	$1.049 \cdot 10^{- 4}$	$2.236 \cdot 10^{- 5}$	$1.528 \cdot 10^{- 5}$	$1.225 \cdot 10^{- 5}$	$1.049 \cdot 10^{- 5}$

To have clear information on the accuracy of the simulated results, in Tables 4 Table 5–6 the deviation of the MC results from the benchmarks of Sec. 2.1 is shown. The deviation is normalized to the standard error of the MC results (this defines the $t$ random variable to be used for the one-sample $t$ -test). We found that all the MC results (with only one exception) were consistent with the RTE values within two standard errors. The $t$ -tests for all the cases reported in Tables 4–6 at a 5% level of significance show that the null hypothesis is not rejected. These results testify to the consistency of the extracted trajectories with the statistical rules of propagation of photon transport. The expected values from the RTE are shown in Table 7. The MC code is also able to describe the small differences between the moments for pure isotropic scattering (HG, $g = 0$ ) and Rayleigh scattering as can be noted in the reported tables.

Table 4

Deviation of the MC results from the theory normalized to the MC standard error, ⟨…⟩MC−⟨…⟩RTEσ⟨…⟩MC, for the statistical moments of the coordinates of scattering events in a non-absorbing infinite medium with scattering coefficient μs=1 mm−1, isotropic scattering function obtained with the HG model (g=0, g2=1/3), and three values N of generated trajectories.

⟨…⟩MC−⟨…⟩RTEσ⟨…⟩MC	N=106				N=108				N=1010
k	1	2	3	4	1	2	3	4	1	2	3	4
$⟨ x_{k} ⟩$		−0.865	−0.483	0.271		−0.837	−0.480	−0.491		0.996	0.838	0.467
$⟨ y_{k} ⟩$		0.461	0.052	0.532		−0.257	−0.096	0.798		1.046	−0.089	−0.198
$⟨ z_{k} ⟩$	0.848	−0.162	−0.380	0.023	−0.817	−0.261	0.275	1.022	1.141	1.687	1.734	1.459
$⟨ x_{k}^{2} ⟩$		−1.050	−0.051	−0.362		−0.817	1.244	0.118		−0.746	−0.328	−0.193
$⟨ y_{k}^{2} ⟩$		−1.165	−0.642	−0.710		−0.023	−0.628	−0.244		0.561	0.538	0.093
$⟨ z_{k}^{2} ⟩$	0.981	0.022	0.081	0.190	−0.878	−0.430	0.410	0.235	1.227	1.093	1.769	1.312
$⟨ ρ_{k}^{2} ⟩$		−1.394	−0.378	−0.692		−0.528	0.394	−0.082		−0.116	0.134	−0.065
$⟨ d_{k}^{2} ⟩$	0.982	−0.667	−0.159	−0.300	−0.878	−0.608	0.531	0.113	1.227	0.827	1.380	0.874
$⟨ l_{k} ⟩$	0.848	−0.933	−0.400	−0.030	−0.817	−1.470	−0.191	−0.384	1.141	0.602	0.053	−0.012
$⟨ l_{k}^{2} ⟩$	0.982	−0.525	−0.169	−0.085	−0.878	−1.410	−0.258	0.043	1.227	0.580	0.219	0.048

Table 5

Deviation of the MC results from the theory normalized to the MC standard error, ⟨…⟩MC−⟨…⟩RTEσ⟨…⟩MC, for the statistical moments of the coordinates of scattering events in a non-absorbing infinite medium with scattering coefficient μs=1 mm−1, an anisotropic scattering function obtained with the HG model with g=0.9 (g2=2.62/3), and three values N of generated trajectories.

⟨…⟩MC−⟨…⟩RTEσ⟨…⟩MC	N=106				N=108				N=1010
k	1	2	3	4	1	2	3	4	1	2	3	4
$⟨ x_{k} ⟩$		−0.438	0.709	0.635		0.632	1.750	1.140		0.184	0.401	0.825
$⟨ y_{k} ⟩$		0.714	−0.066	−0.450		−1.130	−1.869	−1.432		0.218	0.479	0.188
$⟨ z_{k} ⟩$	0.848	−1.037	−0.854	−0.450	−0.817	−1.293	0.188	0.103	1.141	0.543	0.175	0.310
$⟨ x_{k}^{2} ⟩$		0.493	0.158	0.527		−0.965	−1.323	−1.026		−1.717	−1.467	−0.280
$⟨ y_{k}^{2} ⟩$		−0.021	−0.072	−0.240		−0.632	−0.189	−0.029		−0.807	−0.836	−0.932
$⟨ z_{k}^{2} ⟩$	0.982	−0.615	−0.275	−0.162	−0.878	−1.289	−0.050	0.071	1.227	0.693	0.703	0.611
$⟨ ρ_{k}^{2} ⟩$		0.294	0.0538	0.180		−0.985	−0.942	−0.665		−1.556	−1.435	−0.764
$⟨ d_{k}^{2} ⟩$	0.982	−0.558	−0.251	−0.104	−0.878	−1.390	−0.244	−0.110	1.227	0.459	0.371	0.370
$⟨ l_{k} ⟩$	0.848	−0.933	−0.400	−0.030	−0.817	−1.470	−0.191	−0.384	1.141	0.602	0.053	−0.012
$⟨ l_{k}^{2} ⟩$	0.982	−0.526	−0.169	−0.085	−0.878	−1.410	−0.258	0.043	1.227	0.580	0.219	0.048

Table 6

Deviation of the MC results from the theory normalized to the MC standard error, ⟨…⟩MC−⟨…⟩RTEσ⟨…⟩MC, for the statistical moments of the coordinates of scattering events in a non-absorbing infinite medium with scattering coefficient μs=1 mm−1, a Rayleigh scattering function (g=0, g2=0.4), and three values N of generated trajectories.

⟨…⟩MC−⟨…⟩RTEσ⟨…⟩MC	N=106				N=108				N=1010
k	1	2	3	4	1	2	3	4	1	2	3	4
$⟨ x_{k} ⟩$		−0.845	0.623	−0.092		−0.848	−0.964	−0.294		1.108	0.198	0.393
$⟨ y_{k} ⟩$		0.044	−0.035	0.172		−0.235	0.072	−0.069		1.096	0.447	−0.308
$⟨ z_{k} ⟩$	0.848	1.451	0.996	1.122	−0.817	−0.943	−0.382	−0.397	1.141	−0.031	0.268	0.555
$⟨ x_{k}^{2} ⟩$		−1.044	−1.295	−0.723		−0.709	0.542	−0.556		−0.660	−1.455	−0.701
$⟨ y_{k}^{2} ⟩$		−1.158	−0.0546	0.962		0.002	0.063	−0.142		0.486	−1.767	−1.660
$⟨ z_{k}^{2} ⟩$	0.982	0.554	0.769	0.308	−0.878	−1.691	−0.747	−0.697	1.227	0.446	0.993	0.960
$⟨ ρ_{k}^{2} ⟩$		−1.384	−1.170	−1.090		−0.444	0.386	−0.452		−0.109	−2.058	−1.527
$⟨ d_{k}^{2} ⟩$	0.982	−0.164	−0.067	−0.441	−0.878	−1.604	−0.347	−0.769	1.227	0.320	−0.395	−0.245
$⟨ l_{k} ⟩$	0.848	−0.933	−0.400	−0.030	−0.817	−1.470	−0.191	−0.384	1.141	0.602	0.053	−0.012
$⟨ l_{k}^{2} ⟩$	0.982	−0.526	−0.169	−0.086	−0.878	−1.410	−0.258	0.043	1.227	0.580	0.219	0.048

Table 7

Expected values from the RTE theory of the statistical moments of the coordinates of scattering events for a non-absorbing infinite medium with scattering coefficient μs=1 mm−1 and three scattering functions: HG (g=0), HG (g=0.9), and Rayleigh.

RTE	HG (g=0)				HG (g=0.9)				Rayleigh
k	1	2	3	4	1	2	3	4	1	2	3	4
$⟨ x_{k} ⟩$ (mm)	0.	0.	0.	0.	0.	0.	0.	0.	0.	0.	0.	0.
$⟨ y_{k} ⟩$ (mm)	0.	0.	0.	0.	0.	0.	0.	0.	0.	0.	0.	0.
$⟨ z_{k} ⟩$ (mm)	1.	1.	1.	1.	1.	1.9	2.71	3.349	1.	1.	1.	1.
$⟨ x_{k}^{2} ⟩$ ( ${mm}^{2}$ )	0.	$0. \bar{6}$	$1. \bar{3}$	2.	0.	$0.12 \overline{6}$	0.469933	1.091246	0.	0.6	1.26	1.926
$⟨ y_{k}^{2} ⟩$ ( ${mm}^{2}$ )	0.	$0. \bar{6}$	$1. \bar{3}$	2.	0.	$0.12 \overline{6}$	0.469933	1.091246	0.	0.6	1.26	1.926
$⟨ z_{k}^{2} ⟩$ ( ${mm}^{2}$ )	2.	$2. \bar{6}$	$3. \bar{3}$	4.	2.	$5.54 \overline{6}$	10.28013	15.91551	2.	2.8	3.48	4.148
$⟨ ρ_{k}^{2} ⟩$ ( ${mm}^{2}$ )	0.	$1. \bar{3}$	$2. \bar{6}$	4.	0.	$2.25 \overline{3}$	0.939866	2.182492	0.	1.2	2.52	3.852
$⟨ d_{k}^{2} ⟩$ ( ${mm}^{2}$ )	2.	4.	6.	8.	2.	5.8	11.22	18.098	2.	4.	6.	8.
$⟨ l_{k} ⟩$ (mm)	1.	2.	3.	4.	1.	2.	3.	4.	1.	2.	3.	4.
$⟨ l_{k}^{2} ⟩$ ( ${mm}^{2}$ )	2.	6.	12.	20.	2.	6.	12.	20.	2.	6.	12.	20.

For higher scattering orders ( $k > 4$ ), we can still use Eqs. (25)–(31). As an example in Table 8, we show the deviation of the MC results for the moment $⟨ d_{k}^{2} ⟩$ from Eq. (31) up to the tenth order for the different scattering functions considered. Except for one case, the deviation between MC results and reference RTE expected values is within two standard errors. These results confirm the same level of accuracy obtained in the other presented tables with the benchmarks used for $k \in {1, \dots, 4}$ . Similar results have also been obtained for the moments $⟨ x_{k}^{2} ⟩$ , $⟨ y_{k}^{2} ⟩$ , and $⟨ z_{k}^{2} ⟩$ for the case of isotropic scattering [Eqs. (28) and (29)].

Table 8

Deviation of the MC results from the theory normalized to the MC standard error, ⟨dk2⟩MC−⟨dk2⟩RTEσ⟨dk2⟩MC, for the statistical moment ⟨dk2⟩ in a non-absorbing infinite medium with scattering coefficient μs=1 mm−1 for three different scattering functions: HG scattering function with g=0 and g=0.9, and Rayleigh scattering function. The number of simulated trajectories is N=1010.

⟨dk2⟩MC−⟨dk2⟩RTEσ⟨dk2⟩MC	N=1010
k	1	2	3	4	5	6	7	8	9	10
HG $g = 0$	−0.150	−1.306	−0.931	−0.386	−0.785	−0.520	−0.312	−0.351	−0.756	−0.308
HG $g = 0.9$	1.318	1.635	2.650	1.643	1.161	0.231	0.418	−0.017	0.239	0.379
Rayleigh	−0.117	−1.336	−0.687	−1.683	−1.548	−1.058	−1.212	−1.054	−1.057	−0.985

3.2.

Second Step of the Verification (Sec. 2.2)

In this section, the radiance, the fluence rate, and the total mean path length spent inside a layered non-absorbing slab subjected to a Lambertian illumination calculated with MC simulations are compared with the exact analytical solutions of Sec. 2.2. This kind of verification is aimed to test the MC code regarding the effects of boundaries on photon migration. To this purpose, we have considered a layered slab geometry. In the considered examples, the boundary effects are tested both internally to the medium, betwen different layers, and also for the boundary with the external medium.

For this second step of the verification method, we have used Eqs. (33) and (37) for the radiance, Eqs. (34) and (38) for the fluence rate, and Eqs. (35) and (39) for the mean path lengths. We notice that, to the best of our knowledge, the methodical application of the radiance and the fluence rate for a verification procedure of MC codes is original. Radiance and fluence are fundamental quantities of transport theory and their correct calculation is crucial in many applications and also for the verification of MC codes.⁴⁰^,⁴¹

3.2.1.

Comparisons for the radiance

The MC simulations have been carried out assuming a unitary incoming flux. In the solutions given in Sec. 2.2, this is equivalent to assume in the theory $I_{0} = \frac{1}{π}$ $W m^{- 2} {sr}^{- 1}$ . All the layered slabs considered in this section are laterally infinitely extended and subjected to a uniform Lambertian illumination.²⁹ The required uniform isotropic radiance illumination (Lambertian) of the slab can be carried out by using the intrinsic symmetries of this geometry that allow to replace the uniform illumination on the external surface by a single point illumination. This is possible by making use of the reciprocity theorem (or geometry equivalence) that allows to swap source and detector.⁴² Moreover, for the slab to have uniform illumination on both sides it is also needed to switch alternatively the illumination, i.e., the injection inside the slab of subsequent simulated photons, from one side to the other of the slab. As a first test for the verification of the radiance we considered a four-layered slab, where each layer was 2.5 mm thick and where we had two different distributions of the refractive index, defined as “Up” and “Dw.” These two distributions were characterized by an increasing and decreasing trend across the slab, respectively (Fig. 1). The external refractive index is 1 for the profile Up and 2 for the profile Dw. The scattering coefficient was fixed to $μ_{s} = 1 {mm}^{- 1}$ in each layer of the slab and the scattering function was derived from the HG model with $g = 0$ . In Fig. 2, the radiance calculated with the MC code is compared to the solution of Eq. (33) for two numbers of simulated trajectories, i.e., $N = 10^{6}$ and $10^{8}$ . The radiance is shown in each layer. The convergence of the simulated radiance to the true values of the invariant solution can be clearly noted in each layer and for both distributions of the refractive index. In an error-free code, the improvement of the accuracy versus $N$ can be also easily visualized for the radiance, that is isotropic and constant inside any sub-volume with a constant refractive index (see Sec. 2.2). This is shown in Fig. 2 where the fluctuations of the simulated radiance around the true values decrease for the larger $N$ . In Fig. 2, it can also be noted the larger oscillations around 90 deg caused by the reduced sampling of photons around this angle.

Fig. 1

Profiles of the refractive index Up and Dw inside a four-layered slab of thickness 10 mm used in the MC simulations are shown in Fig. 2. The external refractive index is 1 for the profile Up and 2 for the profile Dw.

Fig. 2

Radiance inside each layer of a non-absorbing four-layered slab of thickness 10 mm. Inside the slab the scattering coefficient is $μ_{s} = 1 {mm}^{- 1}$ and the phase scattering function is obtained from the HG model with $g = 0$ . Two profiles of the refractive index inside the slab have been considered: Up and Dw, see Fig. 1. The external refractive index is 1 for the profile Up and 2 for the profile Dw. The solution of Eq. (33) (red curves) is compared with the results of MC simulations (noisy curves).

We have further tested the behavior of the radiance generated with the MC code by considering 100-layered slab, where each layer was 0.1 mm thick and where we had two different discrete distributions of the refractive index shown in Fig. 3 and denoted as Up, with an increasing refractive index from one side to the other, and Dw, with a decreasing refractive index from one side to the other. Also in this case the external refractive index is 1 for the profile Up and 2 for the profile Dw. For this 100-layered slab we have considered two values of the scattering coefficient: $μ_{s} = 0$ and $1 {mm}^{- 1}$ . The scattering function was the HG model with $g = 0$ . Figures. 4 and 5 show the comparisons of the radiance calculated with MC simulations and with Eqs. (33) ( $μ_{s} \neq 0$ ) and (37) ( $μ_{s} = 0$ ) for the profiles Up and Dw, respectively. The results are shown for four selected depths $z$ inside the slab, i.e., four selected layers, shown in the inset of the figures. For both refractive index profiles and scattering values, the MC simulations match the exact values of the RTE solutions. For the profile Up and $μ_{s} = 0 {mm}^{- 1}$ , we have the conditions of guided propagation and Eq. (33), i.e., $I_{j} = \frac{1}{π} {(\frac{n_{j}}{n_{e}})}^{2} W m^{- 2} {sr}^{- 1}$ , for the radiance must be replaced by Eq. (37). In this case, the discontinuity of the radiance can be visualized for two values of the angle $θ$ , as it is visible in Fig. 4. For incident angles greater than the maximum entrance angle for a given layer, the radiance drops to zero. This behavior is accurately reproduced by the MC results and is observed for all the depths considered. We can thus conclude that the calculations of radiance are widely verified with the exact solutions given in Sec. 2.2.

Fig. 3

Profile of the refractive index Up and Dw inside 100-layered slab of 10 mm thickness used in the MC simulations shown in Figs. 4–9. The external refractive index is 1 for the profile Up and 2 for the profile Dw.

Fig. 4

Radiance inside a non-absorbing 100-layered slab of thickness 10 mm calculated at different depths from the external boundary and for both the non-scattering (a) and scattering case (b). Inside the slab, we considered an HG phase scattering function with $g = 0$ . The profile of refractive index Up is considered, see Fig. 3. The solutions of Eq. (33) [red lines in (b)] and Eq. (37) [red lines in (a)] are compared with MC results (symbols).

Fig. 5

Radiance inside a non-absorbing 100-layered slab of thickness 10 mm calculated at different depths from the external boundary and for both the non-scattering (a) and scattering case (b). Inside the slab, we considered an HG phase scattering function with $g = 0$ . The profile of refractive index Dw is considered, see Fig. 3. The solutions of Eq. (33) [red lines in (b)] and Eq. (37) [red lines in (a)] are compared with MC results (symbols).

3.2.2.

Comparisons for the fluence rate

In Figs. 6 and 7, for the same 100-layered slab used in the previous figures, we extended the verification to the fluence rate inside the slab for several values of the scattering coefficient and with a HG scattering function with $g = 0$ . The agreement between MC results and Eqs. (34) and (38) is excellent for all the values of the scattering coefficient and for both the refractive index distributions Up and Dw. The figures also show (right panels) the deviation between MC results and the RTE solutions which is within about two standard errors of the calculated fluence. The MC results have been reported for a subset of 100 layers to assure a clear reading of the symbols versus the depth $z$ inside the slab. For the profile Up there is a discontinuity of the fluence rate between the scattering and non-scattering cases: the fluence switches from the value given by the invariance property [Eq. (34)], i.e., $Φ_{j} = 4 (\frac{n_{j}}{n_{e}})^{2} W m^{- 2}$ , to the value obtained with Eq. (38). In Fig. 7, the two RTE curves for $μ_{s} > 0$ and $μ_{s} = 0$ are indistinguishable. We can conclude that also the calculations of fluence rate are very well verified with the exact solutions given in Sec. 2.2.

Fig. 6

Fluence rate for the profiles Up inside 100-layered slab for several values of the scattering coefficient. (a) The results of MC simulations (symbols) and the RTE solution (continuous and broken lines). The fluence rate is plotted against the depth from the external surface of the slab. (b) The deviation between simulations and RTE solution normalized to the standard error.

Fig. 7

Fluence rate for the profiles Dw inside 100-layered slab for several values of the scattering coefficient. (a) The results of MC simulations (symbols) and the RTE solution (continuous and broken lines). The fluence rate is plotted against the depth from the external surface of the slab. (b) The deviation between the Monte Carlo simulations and the RTE solution normalized to the standard error.

3.2.3.

Comparisons for the mean path length

In Figs. 8 and 9, for the same 100-layered slab used for generating Figs. 4–7, we compared the partial mean path length $⟨ L_{j} ⟩$ spent inside the layers of the slab calculated with MC simulations and with Eqs. (35) and (39). The agreement between MC and exact solutions is excellent and very similar to that obtained for the fluence rate in the previous figures. For the profile Up, we also have a discontinuity of $⟨ L_{j} ⟩$ between $μ_{s} \neq 0$ and $μ_{s} = 0$ : $⟨ L_{j} ⟩$ switches from the value given by the invariance property [Eq. (35)], i.e., $⟨ L_{j} ⟩ = 2 s_{j} (\frac{n_{j}}{n_{e}})^{2}$ ( $s_{j}$ thickness of the layer $j$ , in this case, $s_{j} = 0.1 mm$ ), to the value obtained with Eq. (39). In Fig. 9, the two RTE curves for $μ_{s} > 0$ and $μ_{s} = 0$ are indistinguishable. The deviation between MC results and RTE solutions is also in this case within about two standard errors. We have thus evidence that also the calculations of the mean path length are widely verified with the exact solutions given in Sec. 2.2.

Fig. 8

Average internal path length $⟨ L_{j} ⟩$ spent in the layers of 100-layered slab for the profiles Up for several values of the scattering coefficient. (a) The results of MC simulations (symbols) and the RTE solution (continuous and broken lines). The path length is plotted against the depth from the external surface of the slab. (b) The deviation between the Monte Carlo simulations and the RTE solution normalized to the standard error.

Fig. 9

Average internal path length $⟨ L_{j} ⟩$ spent in the layers of 100-layered slab for the profiles Dw for several values of the scattering coefficient. (a) The results of MC simulations (symbols) and the RTE solution (continuous and broken lines). The path length is plotted against the depth from the external surface of the slab. (b) The deviation between simulations and RTE solution normalized to the standard error.

4. Discussion

In the proposed two-step verification method, we have exploited the complementary characteristics of two kinds of benchmarks: one set of analytical benchmarks sensitive to the single-scattering properties in an infinite medium and another set sensitive to the effects of boundaries. Therefore, the two steps are useful to test an MC code for the correctness of trajectories extraction (first step) and the correctness of photons intersection with boundaries, including the calculations of partial and total path lengths (second step). The proposed method significantly extends the verification methods previously published in Refs. 2 and 3. In fact, the case of isotropic scattering in an infinite medium can now be applied for any scattering order [Eqs. (28)–(31)]. Moreover, the benchmark proposed in Martelli et al.,³ which specifically focused on the calculation of partial and total mean path lengths, now also includes other fundamental radiometric quantities such as radiance and fluence rate.

The presented results show an increasing accuracy with the number of simulated trajectories, which is limited only by machine precision. In practice, the results and the related $t$ -tests confirm the reliability of our MC code and serve as an example of the proposed approach.

Probably, at this point, it is not useless to repeat the following concept: the success of a test (i.e., a null hypothesis not rejected) does not guarantee that a particular MC code performs correctly for all the cases envisioned. Only the failure of the test gives the certainty of having an error inside the code.¹ Thus, the careful reader should realize that, from a practical point of view, no existing method allows us to detect 100% of the errors in an MC code. Nevertheless, we believe that our proposed method is a valuable contribution to increasing the sensitivity of a test to detect coding errors. Also, the two-step method allows one to be more specific about the origin of the error, i.e., if it is caused by an incorrect implementation of the phase function and trajectories’ extraction (step 1) or by an incorrect treatment of the boundaries (step 2). Thus, in this frame, the present contribution certainly represents a fundamental improvement.

Another observation is about testing for absorption effects. Biological tissues are absorbing media and the absorption coefficient is always considered in MC simulations. We notice that even if the two-step method involves only non-absorbing media, this fact does not compromise the generality of the verification. In fact, absorption does not change the photons’ trajectories and its effect can be accounted for with a straightforward implementation of the mBLL,²⁸ without any computational criticality. At the core of the application of the mBLL is the correct evaluation of the partial path lengths, which is tested in step 2 of our method.

One final observation is about the extension of step 2. In the examples reported in this work, we simplified the application of the second step by using the symmetries of the media considered, which allowed us to implement the Lambertian illumination at one point of the external boundary and use the reciprocity theorem to calculate the quantities of interest. However, we stress that the proposed method can be applied to any complex geometry where a full Lambertian illumination can be implemented. This would allow one to test an MC code (for the boundary effects) directly for more realistic cases of interest, without resorting to simplifying the assumptions on the geometry and the distribution of the scattering properties. Future research will be directed to verify the feasibility of this point.

5. Appendix: Moments for Isotropic Scattering Functions

In this appendix, we prove the validity of Eqs. (28)–(31) in Sec. 2.1.2. Equation (31) can be derived with an exact procedure from the RTE according to the results of Ref. 32. A separate proof must be given for Eqs. (28)–(30). We observe that the distribution function of the random variable $z_{k}^{2}$ , due to the pencil beam source emitting along $z$ and to the hypothesis of isotropic scattering, differs from the distributions of $x_{k}^{2}$ and $y_{k}^{2}$ by the first scattering order. After the first scattering, due to the isotropic distribution of the scattering angle, the two distributions are indistinguishable. For the first scattering order, both distributions of $x_{k}^{2}$ and $y_{k}^{2}$ are null. This fact implies that the difference between $⟨ z_{k}^{2} ⟩$ and $⟨ x_{k}^{2} ⟩$ is given by the right term of Eq. (4). Let’s verify this property. The random variable $z_{k}$ can be written as

Eq. (40)

z_{k} = z_{1} + (z_{k} - z_{1}) = z_{1} + Δ z_{k} .

Therefore,

Eq. (41)

z_{k}^{2} = z_{1}^{2} + Δ z_{k}^{2} + 2 z_{1} Δ z_{k} .

We note that the random variables

z_{1}

and

Δ z_{k}

are independent and also that

⟨ Δ z_{k} ⟩ = 0

because of the hypothesis of isotropic scattering. Thus, taking the average of the above equation, we get

Eq. (42)

⟨ z_{k}^{2} ⟩ = ⟨ z_{1}^{2} ⟩ + ⟨ Δ z_{k}^{2} ⟩ + 2 ⟨ z_{1} ⟩ ⟨ Δ z_{k} ⟩ = ⟨ z_{1}^{2} ⟩ + ⟨ Δ z_{k}^{2} ⟩ .

And, because of the isotropic scattering, the random variables

Δ z_{k}

,

x_{k}

, and

y_{k}

have the same distribution function. Thus we have

Eq. (43)

⟨ z_{k}^{2} ⟩ = ⟨ x_{k}^{2} ⟩ + \frac{2}{μ_{s}^{2}} = ⟨ y_{k}^{2} ⟩ + \frac{2}{μ_{s}^{2}} .

By exploiting Eqs. (31) and (43), we can write

Eq. (44)

3 ⟨ x_{k}^{2} ⟩ + \frac{2}{μ_{s}^{2}} = \frac{2 k}{μ_{s}^{2}},

from which we obtain

Eq. (45)

⟨ x_{k}^{2} ⟩ = ⟨ y_{k}^{2} ⟩ = \frac{2}{3 μ_{s}^{2}} (k - 1),

which proves Eq. (28), and

Eq. (46)

⟨ z_{k}^{2} ⟩ = \frac{2}{3 μ_{s}^{2}} (k + 2),

which proves Eq. (29). Finally, Eq. (30) is a trivial consequence of the above relations. Thus, all the relations given in Sec. 2.1.2 for isotropic scattering are exact within the RTE.

Disclosures

The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.

Acknowledgments

The authors wish to thank Prof. Giovanni Zaccanti for his precious suggestions. The authors wish to acknowledge support from Fondazione Cassa di Risparmio di Firenze, Grant No. 2020, and also the National Institutes of Health, Grants Nos. R01 EB029414 and R01 NS095334.

References

1.

B. D. Ganapol, Analytical Benchmarks for Nuclear Engineering Applications, Nuclear Energy Agency(2008). Google Scholar

2.

G. Zaccanti et al., “Analytic relationships for the statistical moments of scattering point coordinates for photon migration in a scattering medium,” Pure Appl. Opt., 3 897 –905 (1994). https://doi.org/10.1088/0963-9659/3/5/019 PAOAE3 0963-9659 Google Scholar

3.

F. Martelli et al., “Verification method of Monte Carlo codes for transport processes with arbitrary accuracy,” Sci. Rep., 11 19486 (2021). https://doi.org/10.1038/s41598-021-98429-3 SRCEC3 2045-2322 Google Scholar

4.

E. Berrocal et al., “Laser light scattering in turbid media Part I: Experimental and simulated results for the spatial intensity distribution,” Opt. Express, 15 10649 –10665 (2007). https://doi.org/10.1364/OE.15.010649 OPEXFF 1094-4087 Google Scholar

5.

D. Frantz, J. Jönsson and E. Berrocal, “Multi-scattering software Part II: Experimental validation for the light intensity distribution,” Opt. Express, 30 1261 –1279 (2022). https://doi.org/10.1364/OE.445394 OPEXFF 1094-4087 Google Scholar

6.

D. A. Boas et al., “Three-dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express, 10 159 –170 (2002). https://doi.org/10.1364/OE.10.000159 OPEXFF 1094-4087 Google Scholar

7.

Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Biomed. Opt. Express, 17 20179 –20190 (2009). https://doi.org/10.1364/OE.17.020178 BOEICL 2156-7085 Google Scholar

8.

Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express, 1 165 –175 (2010). https://doi.org/10.1364/BOE.1.000165 BOEICL 2156-7085 Google Scholar

9.

H. Shen and G. Wang, “A tetrahedron-based inhomogeneous Monte Carlo optical simulator,” Phys. Med. Biol., 55 947 –962 (2010). https://doi.org/10.1088/0031-9155/55/4/003 PHMBA7 0031-9155 Google Scholar

10.

N. Ren et al., “GPU-based Monte Carlo simulation for light propagation in complex heterogeneous tissues,” Opt. Express, 18 6811 –6823 (2010). https://doi.org/10.1364/OE.18.006811 OPEXFF 1094-4087 Google Scholar

11.

H. Shen and G. Wang, “A study on tetrahedron-based inhomogeneous Monte Carlo optical simulation,” Biomed. Opt. Express, 2 44 –57 (2011). https://doi.org/10.1364/BOE.2.000044 BOEICL 2156-7085 Google Scholar

12.

A. Doronin and I. Meglinski, “Online object oriented Monte Carlo computational tool for the needs of biomedical optics,” Biomed. Opt. Express, 2 2461 –2469 (2011). https://doi.org/10.1364/BOE.2.002461 BOEICL 2156-7085 Google Scholar

13.

R. Watté et al., “Modeling the propagation of light in realistic tissue structures with MMC-fpf: a meshed Monte Carlo method with free phase function,” Opt. Express, 23 17467 –17486 (2015). https://doi.org/10.1364/OE.23.017467 OPEXFF 1094-4087 Google Scholar

14.

J. Cassidy et al., “High-performance, robustly verified Monte Carlo simulation with FullMonte,” J. Biomed. Opt., 23 (8), 085001 (2018). https://doi.org/10.1117/1.JBO.23.8.085001 JBOPFO 1083-3668 Google Scholar

15.

C. J. Zoller et al., “Parallelized Monte Carlo software to efficiently simulate the light propagation in arbitrarily shaped objects and aligned scattering media,” J. Biomed. Opt., 23 065004 (2018). https://doi.org/10.1117/1.JBO.23.6.065004 JBOPFO 1083-3668 Google Scholar

16.

A. Leino, A. Pulkkinen and T. Tarvainen, “ValoMC: a Monte Carlo software and MATLAB toolbox for simulating light transport in biological tissue,” OSA Continuum, 2 957 –972 (2019). https://doi.org/10.1364/OSAC.2.000957 Google Scholar

17.

S. Yan and Q. Fang, “Hybrid mesh and voxel based Monte Carlo algorithm for accurate and efficient photon transport modeling in complex bio-tissues,” Biomed. Opt. Express, 11 6262 –6270 (2020). https://doi.org/10.1364/BOE.409468 BOEICL 2156-7085 Google Scholar

18.

L. Wang, S. Ren and X. Chen, “Comparative evaluations of the Monte Carlo-based light propagation simulation packages for optical imaging,” J. Innov. Opt. Health Sci., 11 1750017 (2018). https://doi.org/10.1142/S1793545817500171 Google Scholar

19.

L. Wang, S. L. Jacques and L. Zheng, “MCML Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Prog. Biomed., 47 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F Google Scholar

20.

A. Liemert and A. Kienle, “Analytical Green’s function of the radiative transfer radiance for the infinite medium,” Phys. Rev. E, 83 036605 (2011). https://doi.org/10.1103/PhysRevE.83.036605 Google Scholar

21.

A. Liemert and A. Kienle, “Light transport in three-dimensional semi-infinite scattering media,” J. Opt. Soc. Am. A, 29 1475 –1481 (2012). https://doi.org/10.1364/JOSAA.29.001475 JOAOD6 0740-3232 Google Scholar

22.

A. Liemert and A. Kienle, “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep., 3 2018 (2013). https://doi.org/10.1038/srep02018 SRCEC3 2045-2322 Google Scholar

23.

A. Liemert, D. Reitzle and A. Kienle, “Analytical solutions of the radiative transport equation for turbid and fluorescent layered media,” Sci. Rep., 7 3819 (2017). https://doi.org/10.1038/s41598-017-02979-4 SRCEC3 2045-2322 Google Scholar

24.

H. C. van de Hulst, Light Scattering by Small Particles, John Wiley and Sons, New York (1957). Google Scholar

25.

H. C. van de Hulst, Multiple Light Scattering: Tables, Formulas, and Applications, Academic Press, New York (1980). Google Scholar

26.

R. G. Giovanelli, “Reflection by semi-infinite diffusers,” Opt. Acta, 2 153 –162 (1955). https://doi.org/10.1080/713821040 OPACAT 0030-3909 Google Scholar

27.

L. V. Wang and H. Wu, Biomedical Optics, Principles and Imaging, John Wiley and Sons, New York (2007). Google Scholar

28.

F. Martelli et al., Light Propagation through Biological Tissue and Other Diffusive Media: Theory, Solutions, and Software, SPIE Press, Bellingham, Washington (2009). Google Scholar

29.

F. Martelli et al., “Invariance properties of exact solutions of the radiative transfer equation,” J. Quantum Specstrosc. Radiat. Transf., 276 107887 (2021). https://doi.org/10.1016/j.jqsrt.2021.107887 Google Scholar

30.

F. Tommasi et al., “Invariance property in inhomogeneous scattering media with refractive index mismatch,” Phys. Rev. A, 102 043501 (2020). https://doi.org/10.1103/PhysRevA.102.043501 Google Scholar

31.

F. Tommasi et al., “Invariance property in scattering media and absorption,” Opt. Commun., 458 124786 (2020). https://doi.org/10.1016/j.optcom.2019.124786 OPCOB8 0030-4018 Google Scholar

32.

A. Liemert, F. Martelli and A. Kienle, “Analytical solutions for the mean square displacement derived from transport theory,” (2022). Google Scholar

33.

S. Blanco and R. Fournier, “An invariance property of diffusive random walks,” Europhys. Lett., 61 168 –173 (2003). https://doi.org/10.1209/epl/i2003-00208-x EULEEJ 0295-5075 Google Scholar

34.

A. Mazzolo, C. de Mulatier and A. Zoia, “Cauchy’s formulas for random walks in bounded domains,” J. Math. Phys., 55 (8), 083308 (2014). https://doi.org/10.1063/1.4891299 Google Scholar

35.

M. Majic, W. R. C. Somerville and E. C. Le-Ru, “Mean path length inside nonscattering refractive objects,” Phys. Rev. A, 103 L031502 (2021). https://doi.org/10.1103/PhysRevA.103.L031502 Google Scholar

36.

P. Bruscaglioni et al., “Modified Monte Carlo method to evaluate multiple scattering effects on light beam transmission through a turbid atmosphere,” Proc. SPIE, 0369 164 –173 (1983). https://doi.org/10.1117/12.934363 PSISDG 0277-786X Google Scholar

37.

G. Battistelli et al., “Use of two scaling relations in the study of multiple scattering effect on the transmittance of light beams through a turbid atmosphere,” J. Opt. Soc. Am. A, 2 903 –912 (1985). https://doi.org/10.1364/JOSAA.2.000903 JOAOD6 0740-3232 Google Scholar

38.

G. Zaccanti, “Monte Carlo study of light propagation in optically thick media: Point-source case,” Appl. Opt., 30 2031 –2041 (1991). https://doi.org/10.1364/AO.30.002031 APOPAI 0003-6935 Google Scholar

39.

F. Martelli et al., “Photon migration through a turbid slab described by a model based on diffusion approximation. II. Comparison with Monte Carlo results,” Appl. Opt., 36 4600 –4612 (1997). https://doi.org/10.1364/AO.36.004600 APOPAI 0003-6935 Google Scholar

40.

A. Liemert and A. Kienle, “Analytical solution of the radiative transfer equation for infinite-space fluence,” Phys. Rev. A, 83 015804 (2011). https://doi.org/10.1103/PhysRevA.83.015804 Google Scholar

41.

A. Liemert and A. Kienle, “Green’s function of the time-dependent radiative transport equation in terms of rotated spherical harmonics,” Phys. Rev. E, 86 036603 (2012). https://doi.org/10.1103/PhysRevE.86.036603 Google Scholar

42.

A. F. Bielajew, “Fundamentals of the Monte Carlo method for neutral and charged particle transport,” (1998) http://websites.umich.edu/nersb590/CourseLibrary/MCbook.pdf Google Scholar

Biography

Angelo Sassaroli graduated from the University of Florence, Italy, in 1996 (laureate in physics). He received his PhD in physics from the University of Electro-Communications, Tokyo, Japan, in 2002. From July 2002 to August 2007, he was a research associate in the group of Prof. Sergio Fantini at Tufts University, Medford (MA). In September 2007, he was appointed by Tufts University as a research assistant professor. His fields of research are near-infrared spectroscopy and diffuse optical tomography.

Federico Tommasi graduated from the University of Florence, Italy, in 2011. He received his PhD in physics from the Department of Physics and Astronomy of the University of Florence, Italy, in 2015. He is currently in the Department of Physics and Astronomy of the University of Florence. His major topics are random laser sources and their applications together with studies of light propagation through complex active and passive random media by means of numerical simulations.

Stefano Cavalieri received his PhD in physics from the University of Florence, Italy, in 1990. He is appointed as an associate professor in physics matter in the Department of Physics and Astronomy of the University of Florence. His major topics are random laser sources, random media, coherent control of laser–matter interaction, stochastic phenomena, and their applications.

Lorenzo Fini is appointed as a researcher in the Department of Physics and Astronomy of the University of Florence, Italy. His major topic is applied optics.

André Liemert received his Dr.-Ing. degree in electrical engineering from the University of Ulm in 2011. Since 2011, he has worked as a postdoctoral researcher in the Materials Optics and Imaging Group of Prof. Alwin Kienle at the Institut für Lasertechnologien in der Medizin und Meßtechnik an der Universität Ulm.

Alwin Kienle is the speaker of the board of directors at ILM (Institut für Lasertechnologien in der Medizin und Meßtechnik an der Universität Ulm), Ulm, Germany, and head of the department Quantitative Imaging and Sensors at ILM. In addition, he is a professor at the University of Ulm. He studied physics and received his doctoral and habilitation degrees from the University of Ulm. As a postdoc, he worked with research groups in Hamilton, Canada, and in Lausanne, Switzerland.

Tiziano Binzoni: Biography is not available.

Fabrizio Martelli graduated from the University of Florence, Italy, in 1996. He received his PhD in physics from the University of Electro-Communications, Tokyo, Japan, in 2002. He is currently appointed as researcher in the Department of Physics and Astronomy of the University of Florence, Italy. His major topic is the modeling of light propagation through highly scattering media based on the radiative transfer equation and on the diffusion equation.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Angelo Sassaroli, Federico Tommasi, Stefano Cavalieri, Lorenzo Fini, André Liemert, Alwin Kienle, Tiziano Binzoni, and Fabrizio Martelli "Two-step verification method for Monte Carlo codes in biomedical optics applications," Journal of Biomedical Optics 27(8), 083018 (20 April 2022). https://doi.org/10.1117/1.JBO.27.8.083018

Received: 30 December 2021; Accepted: 21 March 2022; Published: 20 April 2022

Access the abstract

JOURNAL ARTICLE
25 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY