An Integral-Equation-Based Variance Reduction Method for Accelerated Monte Carlo Simulations

Abstract

In this work, we introduce a novel variance reduction approach utilising the integral formulation of the radiative transfer equation to calculate the radiance in a planar symmetric slab geometry. Due to its integral nature, our method offers a fundamental advantage over well-established variance reduction methods such as the local estimate technique. As opposed to the local estimate procedure, photons add to the overall radiance not only at specific points of interaction but also throughout each consecutive path element; hence, our variance reduction approach can be thought of as an integral local estimate method. This facilitates a substantial enhancement in statistical efficiency, especially in scenarios where only a small number of scattering events or a high attenuation along the detection paths is to be anticipated. To evaluate the overall performance of the integral approach, we incorporated it into a self-developed GPU-accelerated Monte Carlo software, together with a conventional local estimate implementation adapted to slab geometry for a comprehensive comparison.

1. Introduction

The fraction of radiant power reflected or transmitted by a sample in a specific direction is given by the angular reflectance (AR) or transmittance (AT). Consequently, these quantities are relevant in a wide range of research areas where the interaction of light with scattering media is of interest, including atmospheric science [1,2,3,4], radiation therapy [5,6,7], photodynamic therapy [8], forest surface or soil studies [9,10,11], and computer graphics [12,13], in the form of the bidirectional reflectance distribution function (BRDF). In this context, the Monte Carlo (MC) method, which provides a numerical solution to the radiative transfer equation (RTE), is often the tool of choice for theoretical investigations, as it is considered the gold standard for numerical modelling of light propagation in turbid media.

In many cases, when the MC approach is used to simulate the AR or AT, the signal is not detected in discrete directions, as ideally desired, but rather in extended angular bins [3,4,11,14,15,16]. This implies that the signal is always obtained as a sum over an angular range and deviations from discrete detection angles are to be anticipated. This can be partially compensated for by a higher angular resolution, but the variance will increase as hardly any energy packets, called photons in MC terminology, contributes to the individual detection directions. To this end, a number of variance reduction methods (VRMs) [17,18,19] have been introduced by the Monte Carlo community to address these types of problems. A well-established VRM technique is the local estimate (LE) method, also called last flight estimation, the method of expected values and the method of statistical evaluation, which, by virtue of its simplicity and intuitiveness, is typically used to improve the statistics when the probability of a photon contributing to the signal is close to zero due to the detector’s limited spatial extent, or when the signal is detected within a limited field of view. The LE method involves computing the probability of a photon being scattered in the direction of the detector at each point of interaction. Thus, each scattering event adds a contribution to the overall signal. However, in cases of thin layers or low scattering media, only a limited number of LE contributions are obtained due to the small set of interaction points. Additionally, it is critical to consider the attenuation along the detection path when calculating LE weights. For large angles with respect to the surface normal, this results in inadequate statistics for individual scattering events due to the high level of attenuation along the extended detection paths.

In this work, we present a novel type of VRM, derived from the integral form of the RTE, which, similar to the LE method, allows the calculation of the contributions for arbitrary detection directions, i.e., the contributions to all angles of the AT or AR, during the propagation of a photon. In this case, however, the calculation of the radiance contributions extends beyond discrete interaction points to integrally account for each individual path section, which is why our VRM can also be regarded as an integral LE (iLE) method. The incorporation of entire path segments thus allows a significant variance improvement compared to the LE method, especially in the cases specified above. The sole drawback of this type of approach is that it is limited to geometries for which analytical expressions can be derived, namely, in the case presented here, the planar symmetric radiance calculation for a slab of infinite lateral extent, as, while there are several analytical solutions available for the RTE to calculate the planar radiance, for cases of extreme optical properties, a high computational cost must often be incurred to ensure the accuracy and stability of results. To assess the efficacy of our novel approach, we integrated the iLE method into a self-developed GPU (graphics processing unit)-accelerated MC software. Additionally, we adapted the LE technique to slab geometry and incorporated it into the same MC framework for a thorough comparison.

The rest of the paper is structured as follows. In Section 2, the theoretical framework of the LE and iLE methods is introduced together with the algorithmic implementation by means of pseudo-code. In Section 3, the results are presented and discussed; in particular, the benchmark testing and comparison of the two methods based on selected examples are considered. To conclude, the primary results of this study are summarised in Section 4.

2. Methodology

2.1. Local Estimate Method

To reduce the variance of an MC simulation applying the LE method, at each point of interaction, i.e., each scattering event, the probability that a photon is scattered to a given target position and reaches the detector ballistically without any further interactions is calculated. Accordingly, each photon at each scattering event contributes to the detector signal with an adjusted weight that takes into account the scattering function of the medium and the attenuation along the detection path. For spatially resolved signals, the calculation of the local estimation weights is straightforward and can be found in the literature [19,20]. However, when using the MC method to simulate the radiance at the surface, which in turn can be transferred to the AR or AT, considering a laterally infinite layered geometry as shown in Figure 1, any multiple reflections that could contribute to the signal of a particular detection direction must be included in the calculation of the LE weights in the case of a refractive index mismatch ($n_1\ne n_2$). The composition of the LE weights for this particular case is illustrated in Figure 1.

Consider a photon at location ${\overrightarrow{r}}_1$, with weight $w_1$ and normalised direction of propagation ${\widehat{s}}_p$, which, following the well-known MC methodology [21,22], travels along a free path

(1)$\ell =-{\displaystyle \frac{\ln \zeta }{{\mu }_s}}$

to location ${\overrightarrow{r}}_2$, where $\zeta \in [0,1]$ is a random number drawn from a uniform distribution and ${\mu }_s$ is the scattering coefficient of the medium under consideration. Therefore, the initial weight $w_1$ is reduced to

(2)$w_2=w_1\exp \left(-{\mu }_a\ell \right),$

where ${\mu }_a$ is the absorption coefficient. At this point, the first “direct” LE contribution $w_s$, the upper black solid line, to the radiance in the detection direction $\widehat{s}$ on the transmission side ($z=L$) is

(3)$w_s=w_{2}p\left({\widehat{s}}_p·\widehat{s}\right)\exp \left(-{\mu }_t{\displaystyle \frac{L-z_2}{\left|\mu \right|}}\right)=w_{1}p\left({\widehat{s}}_p·\widehat{s}\right)\exp \left(-{\mu }_t{\displaystyle \frac{L-z_2}{\left|\mu \right|}}-{\mu }_a\ell \right),$

where $\mu =\cos \theta$ and $p\left({\widehat{s}}_p·\widehat{s}\right)$ is the scattering phase function, which gives the probability that a photon propagating in direction ${\widehat{s}}_p$ will be scattered in direction $\widehat{s}$, and ${\mu }_t={\mu }_a+{\mu }_s$ is the extinction coefficient. If we continue to follow the direct path with the weight $w_s$, the blue dashed line, a photon can contribute to the same detection direction again at $z=L$ after having been reflected at the upper boundary layer and subsequently reflected at the lower one, and, therefore, this needs to be included. The corresponding updated weight is $w_s{R}_f^2{e}^2$, where $R_f=R_f\left(\theta \right)$ is the Fresnel reflectance as a function of the polar angle $\theta$ and

(4)$e=e\left(\mu \right)=\exp \left(-{\mu }_t{\displaystyle \frac{L}{\left|\mu \right|}}\right),$

accounts for the attenuation along the reflected path. Similarly, we can calculate the contributions of all subsequent reflections (Figure 1). The sum of all the contributions of the direct path can thus be summarised and calculated via a geometric series

(5)$W_s\left(\mu \right)=w_s\sum _{k=0}^{\infty }{\left(R_f^2{e}^2\right)}^k={\displaystyle \frac{w_s}{1-R_f^2{e}^2}}.$

In addition to the direct path, a photon at the interaction point ${\overrightarrow{r}}_2$ could also hit the lower boundary layer ($z=0$) at the angle $\theta$, the black dashed line, and after (multiple) reflection(s) it may again contribute to the same detection direction at $z=L$. Similarly, the sum of the individual contributions can be determined via a geometric series

(6)${\tilde{W}}_s\left(\mu \right)={\tilde{w}}_s{R}_{f}e\sum _{k=0}^{\infty }{\left(R_f^2{e}^2\right)}^k={\displaystyle \frac{{\tilde{w}}_s{R}_{f}e}{1-R_f^2{e}^2}},$

where

(7)${\tilde{w}}_s=w_{1}p\left({\widehat{s}}_p·\tilde{s}\right)\exp \left(-{\mu }_t{\displaystyle \frac{z_2}{\left|\mu \right|}}-{\mu }_a\ell \right).$

The total contribution of LE at the interaction point ${\overrightarrow{r}}_2$ to the radiance in the transmission direction $\widehat{s}$ is, therefore,

(8)$I_t\left(\mu \right)={\displaystyle \frac{1}{\left|\mu \right|}}\left(W_s\left(\mu \right)+{\tilde{W}}_s\left(\mu \right)\right)={\displaystyle \frac{1}{\left|\mu \right|}}{\displaystyle \frac{w_s+{\tilde{w}}_s{R}_{f}e}{1-R_f^2{e}^2}}.$

Due to the symmetry present here, the LE contribution

(9)$I_r\left(\mu \right)={\displaystyle \frac{1}{\left|\mu \right|}}{\displaystyle \frac{{\tilde{w}}_s+w_s{R}_{f}e}{1-R_f^2{e}^2}},$

in the reflection direction $\tilde{s}$ can also be calculated simultaneously. This is done by swapping the weights $w_s$ and ${\tilde{w}}_s$ in Equations (5) and (6), as is immediately apparent from the contributions shown in Figure 1 at the lower boundary surface. Finally, to obtain the AT or AR from the radiance, one must consider the reflection loss and refraction at the interface, alongside the consequent modification of the solid angle element [22]

(10)$\begin{array}{c}\hfill T\left({\mu }^{\prime }\right)={\mu }^{\prime }{\left({\displaystyle \frac{n_1}{n_2}}\right)}^2\left[1-R_f\left(\sqrt{1-{(n_1/n_2)}^2(1-{{\mu }^{\prime }}^2)}\right)\right]I_t\left(\sqrt{1-{(n_1/n_2)}^2(1-{{\mu }^{\prime }}^2)}\right),\end{array}$

(11)$\begin{array}{c}\hfill R\left({\mu }^{\prime }\right)={\mu }^{\prime }{\left({\displaystyle \frac{n_1}{n_2}}\right)}^2\left[1-R_f\left(\sqrt{1-{(n_1/n_2)}^2(1-{{\mu }^{\prime }}^2)}\right)\right]I_r\left(-\sqrt{1-{(n_1/n_2)}^2(1-{{\mu }^{\prime }}^2)}\right),\end{array}$

where ${\mu }^{\prime }=\cos {\theta }^{\prime }$ is the respective outer angle after refraction relative to the surface normal (see Figure 1).

2.2. Integral Local Estimate Method

For the integral local estimate method, we first require the integral form of the planar symmetric RTE, which is usually written in integro-differential form as

(12)$\left(\mu \frac{\partial }{\partial z}+{\mu }_t\right)I(z,\widehat{s})={\mu }_s{\int }_{4\pi }p(\widehat{s},{\widehat{s}}^{\prime })I(z,{\widehat{s}}^{\prime })\mathrm{d}{\widehat{s}}^{\prime }+S(z,\widehat{s}),$

where $S(z,\widehat{s})$ is the source term. To that end, we first extend the scattering operator to an integral over all the arguments of the radiance $I(z,\widehat{s})$ as

(13)$\left(\mu \frac{\partial }{\partial z}+{\mu }_t\right)I(z,\widehat{s})={\mu }_s{\int }_{4\pi }{\int }_0^{L}p(\widehat{s},{\widehat{s}}^{\prime })I(z^{\prime },{\widehat{s}}^{\prime })\delta (z-z^{\prime })\mathrm{d}{z}^{\prime }\mathrm{d}{\widehat{s}}^{\prime }+S(z,\widehat{s}).$

To obtain an integral equation for $I(z,\widehat{s})$, we seek a kernel function $K(z,\widehat{s},z^{\prime },{\widehat{s}}^{\prime })$ that satisfies

(14)$\left(\mu \frac{\partial }{\partial z}+{\mu }_t\right)K(z,\widehat{s},z^{\prime },{\widehat{s}}^{\prime })=p(\widehat{s},{\widehat{s}}^{\prime })\delta (z-z^{\prime }),$

under the appropriate boundary conditions for an infinite slab with Fresnel specular reflections at the bottom and top surfaces

(15)$\begin{array}{cc}\hfill K(z=0,\mu ,\varphi ,z^{\prime },{\widehat{s}}^{\prime })& =R_f\left(\mu \right)K(z=0,-\mu ,\varphi ,z^{\prime },{\widehat{s}}^{\prime }),\hfill \end{array}$

(16)$\begin{array}{cc}\hfill K(z=L,-\mu ,\varphi ,z^{\prime },{\widehat{s}}^{\prime })& =R_f\left(\mu \right)K(z=L,\mu ,\varphi ,z^{\prime },{\widehat{s}}^{\prime }).\hfill \end{array}$

Additionally, we also need the response of the propagation operator on the left hand side to the source term $S(z,\widehat{s})$, which is the ballistic or unscattered contribution $I_0(z,\widehat{s})$ to the total radiance, fulfilling

(17)$\left(\mu \frac{\partial }{\partial z}+{\mu }_t\right)I_0(z,\widehat{s})=S(z,\widehat{s}).$

Of course, $I_0$ must satisfy the same boundary conditions that were imposed on K. Then, the complete radiance $I(z,\widehat{s})$ can be written as

(18)$I(z,\widehat{s})={\mu }_s{\int }_{4\pi }{\int }_0^{L}I(z^{\prime },{\widehat{s}}^{\prime })K(z,\widehat{s},z^{\prime },{\widehat{s}}^{\prime })\mathrm{d}{z}^{\prime }\mathrm{d}{\widehat{s}}^{\prime }+I_0(z,\widehat{s}),$

which is the integral equation we sought. It is easy to see that if both K and $I_0$ satisfy the required boundary conditions, then I will satisfy them as well. The differential Equation (14) for the kernel function $K(z,\widehat{s},z^{\prime },{\widehat{s}}^{\prime })$ under the Fresnel boundary conditions (15) and (16) can be solved in closed form. The solution is given by

(19)$\begin{array}{c}\hfill K(z,\widehat{s},z^{\prime },{\widehat{s}}^{\prime })=\frac{p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)}{\left|\mu \right|}\exp \left(-{\mu }_t\frac{z-z^{\prime }}{\mu }\right)\Theta \left(\frac{z-z^{\prime }}{\mu }\right)+\frac{R_f\left(\right|\mu \left|\right)}{\left|\mu \right|}\frac{\exp \left(-{\mu }_t\frac{z}{\mu }\right)}{1-R_f^2\left(\right|\mu \left|\right)\exp \left(-2{\mu }_t\frac{L}{\left|\mu \right|}\right)}\\ \hfill \times \left\{\begin{array}{cc}p\left(\tilde{s}·{\widehat{s}}^{\prime }\right)\exp \left(-{\mu }_t\frac{z^{\prime }}{\mu }\right)+R_f\left(\mu \right)p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)\exp \left({\mu }_t\frac{z^{\prime }-2L}{\mu }\right),\hfill & \mu >0\hfill \\ R_f\left(\right|\mu \left|\right)p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)\exp \left({\mu }_t\frac{z^{\prime }+2L}{\mu }\right)+p\left(\tilde{s}·{\widehat{s}}^{\prime }\right)\exp \left({\mu }_t\frac{2L-z^{\prime }}{\mu }\right),\hfill & \mu <0\hfill \end{array}\right.,\end{array}$

and the derivation of the integral kernel can be found in Appendix A. Here, $\tilde{s}$ is the “reflected” direction to $\widehat{s}$ (see Figure 1), with

(20)$\tilde{s}·{\widehat{s}}^{\prime }=-\mu {\mu }^{\prime }+\sqrt{1-{\mu }^2}\sqrt{1-{{\mu }^{\prime }}^2}\cos (\varphi -{\varphi }^{\prime }).$

The solution for $I_0$ depends on the exact source term which is used, i.e., for the commonly assumed collimated source given by $S_c(z,\widehat{s})=\delta \left(z\right)\delta (\widehat{s}-{\widehat{s}}_0)$ with incident direction ${\widehat{s}}_0$, the result is

(21)$\begin{array}{c}\hfill I_0(z,\widehat{s})=\frac{1}{\left|\mu \right|}\frac{1}{1-R_f^2\left(\right|\mu \left|\right)\exp \left(-2{\mu }_t\frac{L}{\left|\mu \right|}\right)}\left[\exp \left(-{\mu }_t\frac{z}{\left|\mu \right|}\right)\delta (\widehat{s}-{\widehat{s}}_0)\right.\\ \hfill \left.+R_f\left(\right|\mu \left|\right)\exp \left(-{\mu }_t\frac{2L-z}{\left|\mu \right|}\right)\delta (\widehat{s}-{\tilde{s}}_0)\right],\end{array}$

which can also be deduced from Equation (19). Again, ${\tilde{s}}_0$ is the reflected direction to ${\widehat{s}}_0$. The integral Equation (18) can be treated in two ways. First, one could attempt to solve it directly, circumventing the complicated approximation methods that stem from the original integro-differential form of the RTE [23]. Second, one could insert a given approximation (that does not solve the equation) on the right hand side, which usually leads to a superior expression for the radiance [24]. Here, we use the radiance $\tilde{I}$ predicted by the Monte Carlo simulation, which is a collection of path segments p, each with direction ${\widehat{s}}_p$ between $z_1$ and $z_2$

(22)$\tilde{I}(z^{\prime },{\widehat{s}}^{\prime })=\sum _p\frac{1}{{\mu }_p}\delta ({\widehat{s}}^{\prime }-{\widehat{s}}_p)\left[\Theta (z^{\prime }-z_1)-\Theta (z^{\prime }-z_2)\right]\exp \left(-{\mu }_a\frac{z^{\prime }-z_1}{{\mu }_p}\right),$

where the indices of the z limits were suppressed. One of these path elements can then be inserted into (18) on the right hand side, dropping the ballistic term $I_0$, which may be re-added later. Evaluating this at $z=0$ for $\mu <0$ and at $z=L$ for $\mu >0$ then leads to the contribution of the path element to the radiance at the top and bottom interfaces, which will be called $I_r\left(\mu \right)$ and $I_t\left(\mu \right)$, respectively. Note that the other halves of $I_r\left(\mu \right)$ and $I_t\left(\mu \right)$ with respect to $\mu$ can be computed via the boundary conditions (15) and (16). Immediately performing the angular integration, we get

(23)$\begin{array}{c}\hfill I_r\left(\mu \right)=\frac{{\mu }_s}{\left|\mu \right|{\mu }_p}{\int }_{z_1}^{z_2}\exp \left(-{\mu }_a\frac{z^{\prime }-z_1}{{\mu }_p}\right)\left[p\left(\widehat{s}·{\widehat{s}}_p\right)\exp \left({\mu }_t\frac{z^{\prime }}{\mu }\right)\right.\\ \hfill +\frac{R_f\left(\right|\mu \left|\right)}{1-R_f^2\left(\right|\mu \left|\right)\exp \left(-2{\mu }_t\frac{L}{\left|\mu \right|}\right)}\left(R_f\left(\right|\mu \left|\right)p\left(\widehat{s}·{\widehat{s}}_p\right)\exp \left({\mu }_t\frac{z^{\prime }+2L}{\mu }\right)\right.\\ \hfill \left.\left.+p\left(\tilde{s}·{\widehat{s}}_p\right)\exp \left({\mu }_t\frac{2L-z^{\prime }}{\mu }\right)\right)\right]\mathrm{d}{z}^{\prime },\end{array}$

and

(24)$\begin{array}{c}\hfill I_t\left(\mu \right)=\frac{{\mu }_s}{\left|\mu \right|{\mu }_p}{\int }_{z_1}^{z_2}\exp \left(-{\mu }_a\frac{z^{\prime }-z_1}{{\mu }_p}\right)\left[p\left(\widehat{s}·{\widehat{s}}_p\right)\exp \left(-{\mu }_t\frac{(L-z^{\prime })}{\mu }\right)\right.\\ \hfill +\frac{R_f\left(\right|\mu \left|\right)}{1-R_f^2\left(\right|\mu \left|\right)\exp \left(-2{\mu }_t\frac{L}{\left|\mu \right|}\right)}\left(R_f\left(\right|\mu \left|\right)p\left(\widehat{s}·{\widehat{s}}_p\right)\exp \left(-{\mu }_t\frac{3L-z^{\prime }}{\mu }\right)\right.\\ \hfill \left.\left.+p\left(\tilde{s}·{\widehat{s}}_p\right)\exp \left(-{\mu }_t\frac{L+z^{\prime }}{\mu }\right)\right)\right]\mathrm{d}{z}^{\prime }.\end{array}$

We then introduce two abbreviations $X\left(\mu \right)$ and $Y\left(\mu \right)$, defined as

(25)$X\left(\mu \right):=\frac{1}{\mu }\frac{\exp \left(\frac{{\mu }_t}{\mu }{z}_1\right)-\exp \left(-{\mu }_a\ell +\frac{{\mu }_t}{\mu }{z}_2\right)}{{\mu }_a-{\mu }_t\frac{{\mu }_p}{\mu }},$

and

(26)$Y\left(\mu \right):=\frac{1}{\mu }\frac{\exp \left(\frac{{\mu }_t}{\mu }(L-z_1)\right)-\exp \left(-{\mu }_a\ell +\frac{{\mu }_t}{\mu }(L-z_2)\right)}{{\mu }_a+{\mu }_t\frac{{\mu }_p}{\mu }},$

where $\ell =(z_2-z_1)/{\mu }_p$ is the path length. After evaluating the remaining integrals in (23) and (24), the final results for $I_r\left(\mu \right)$ with $\mu <0$ and $I_t\left(\mu \right)$ with $\mu >0$ can then be written using these abbreviations as

(27)$\begin{array}{c}\hfill I_r\left(\mu \right)=-{\mu }_s\left[p\left(\widehat{s}·{\widehat{s}}_p\right)X\left(\mu \right)+\frac{R_f\left(\right|\mu \left|\right)\exp \left(-\frac{{\mu }_t}{\left|\mu \right|}L\right)}{1-R_f^2\left(\right|\mu \left|\right)\exp \left(-2\frac{{\mu }_t}{\left|\mu \right|}L\right)}\right.\\ \hfill \left.\phantom{\frac{R_f\left(\right|\mu \left|\right)\exp \left(-\frac{{\mu }_t}{\left|\mu \right|}L\right)}{1-R_f^2\left(\right|\mu \left|\right)}}\times \left(p\left(\widehat{s}·{\widehat{s}}_p\right)R_f\left(\right|\mu \left|\right)\exp \left(-\frac{{\mu }_t}{\left|\mu \right|}L\right)X\left(\mu \right)+p\left(\tilde{s}·{\widehat{s}}_p\right)Y\left(\mu \right)\right)\right],\end{array}$

and

(28)$\begin{array}{c}\hfill I_t\left(\mu \right)=-{\mu }_s\left[p\left(\widehat{s}·{\widehat{s}}_p\right)Y(-\mu )+\frac{R_f\left(\right|\mu \left|\right)\exp \left(-\frac{{\mu }_t}{\left|\mu \right|}L\right)}{1-R_f^2\left(\right|\mu \left|\right)\exp \left(-2\frac{{\mu }_t}{\left|\mu \right|}L\right)}\right.\\ \hfill \left.\phantom{\frac{R_f\left(\right|\mu \left|\right)\exp \left(-\frac{{\mu }_t}{\left|\mu \right|}L\right)}{1-R_f^2\left(\right|\mu \left|\right)}}\times \left(p\left(\widehat{s}·{\widehat{s}}_p\right)R_f\left(\right|\mu \left|\right)\exp \left(-\frac{{\mu }_t}{\left|\mu \right|}L\right)Y(-\mu )+p\left(\tilde{s}·{\widehat{s}}_p\right)X(-\mu )\right)\right].\end{array}$

Note that if the observation angles are always taken relative to the outward surface normal, the values of X and Y in (27) and (28) only differ by a sign. However, if the denominators in $X\left(\mu \right)$ or $Y\left(\mu \right)$ given by

(29)$\begin{array}{cc}\hfill d_x& :={\mu }_a-{\mu }_t\frac{{\mu }_p}{\mu }\hfill \end{array}$

(30)$\begin{array}{cc}\hfill d_y& :={\mu }_a+{\mu }_t\frac{{\mu }_p}{\mu }\hfill \end{array}$

are close to zero, large cancellation errors can occur in the numerical evaluation of (25) or (26). In these cases, the expressions for X or Y can be replaced by the expansions

(31)$X\left(\mu \right)=\frac{\exp \left(\frac{{\mu }_t}{\mu }{z}_1\right)}{\mu }\left(\ell -\frac{{\ell }^2{d}_x}{2}+\frac{{\ell }^3{d}_x^2}{6}-O\left(d_x^3\right)\right),$

and

(32)$Y\left(\mu \right)=\frac{\exp \left(\frac{{\mu }_t}{\mu }(L-z_1)\right)}{\mu }\left(\ell -\frac{{\ell }^2{d}_y}{2}+\frac{{\ell }^3{d}_y^2}{6}-O\left(d_y^3\right)\right).$

2.3. Monte Carlo Implementation

For the following benchmark test of the LE and iLE variance reduction methods described in Section 2.1 and Section 2.2, both were integrated into a self-developed GPU-accelerated MC program. The basic simulation routine will now be illustrated in more detail using the pseudo-code provided in Algorithms 1–3. It should be noted that the main purpose of this section is to show the basic algorithmic integration of the variance reduction methods. For the exact technical implementation, the authors refer readers to the Git Repository (https://github.com/dreitzle/radMC/tree/master, accessed on 20 December 2023), where the complete source code can be found. The GPU code was implemented in OpenCL and the host code in both Python and C++.

Starting with the Monte Carlo framework, Algorithm 1 outlines the propagation of a single photon in the geometry depicted in Figure 1. In lines 1–5, the starting position and direction of the photon are set. In the following, it is always assumed that the started photons hit the lower boundary ($z=0$) at an angle of incidence ${\theta }_0$, where ${\theta }_0^{\prime }$ is the angle after the initial refraction at the surface (see Figure 1). Note that only the z-component of the photon position is stored, as this is the only spatial component needed for further calculations. Furthermore, the direction of propagation is described by the cosine and sine terms of the spherical coordinates, as these are needed to calculate the LE and iLE contributions, respectively, thus avoiding a previous coordinate transformation. Prior to the actual propagation loop, the initial photon weight is adjusted in line 6 according to the Fresnel reflectivity. A photon eventually propagates through the medium until its weight falls below a critical value $ϵ$. Russian roulette then decides with a given survival probability q whether the photon with the adjusted weight $w/q$ will continue to propagate or will be terminated with a probability of $1-q$. For the propagation itself, an MC variant was chosen, in which photons are always scattered at interaction points, and the weight is then reduced according to the absorption properties of the medium (compare Section 2.1). Similarly, when interacting with interfaces, photons are always reflected and the weight is adjusted according to the Fresnel equations. It is worth noting that when a reflection event takes place in lines 12–16 or 17–21, there is no immediate change in direction or weight and is initially only monitored using the Boolean variable $reflected$. The reason for this is that the calculation of the radiance contributions in line 22 can then be called at exactly one point in the code, which allows for better synchronisation of the workgroups, achieving better performance. The algorithmic procedure to calculate the radiance contributions using LE and iLE can be seen in Algorithms 2 and 3, respectively.

The LE implementation first examines whether a photon has already gone through one of the boundary layers. Since the method does not include direct detector contributions, there is no calculation of the radiance contributions in this situation. If the interaction point is located in the medium, the reflectivity is computed for a given detection angle ${\theta }_i$ and the weight taking into account the attenuation along the reflection path, as seen in Equation (4). Then the calculation of the contributions is carried out as described in Section 2.1 for each angle ${\theta }_i$ and the related angles ${\varphi }_j$. Continuing with the iLE procedure (Algorithm 3), a similar scheme is adopted; however, it is redundant to test whether the interaction point is outside the medium at the beginning, as the weights are calculated integrally up to the interfaces. Note that inside the second for loop in lines 10–19 both denominators, $d_x$ and $d_y$, are checked for proximity to zero to prevent cancellation errors while calculating X and Y, as explained in Section 2.2. If one of the denominators falls below its respective threshold ${ϵ}_1$ or ${ϵ}_2$, the corresponding X or Y is replaced by its expansion, where ${ϵ}_1={ϵ}_2={10}^{-4}$ has been found to give numerically stable results. Finally, for improved readability, the variables $r_{\mathrm{fac}}$, $c_1$, and $c_2$ are introduced in lines 20–22, from which the respective contributions to the radiance are calculated.

Algorithm 1 Single photon propagation

1:. $z\leftarrow 0$   2:. $\cos \theta \leftarrow \cos {\theta }_0^{\prime }$                           ▹ Direction after initial refraction   3:. $\sin \theta \leftarrow \sqrt{1-{\cos }^2\theta }$   4:. $\cos \varphi \leftarrow \cos {\varphi }_0$   5:. $\sin \varphi \leftarrow \sin {\varphi }_0$   6:. w ← 1 − CalcReflec$\left({\theta }_0^{\prime }\right)$                ▹ Reduced weight according to Fresnel equations   7:. $reflected\leftarrow \mathbf{false}$   8:. while true do   9:.   $\ell \leftarrow -\frac{\ln \zeta }{{\mu }_s}$ 10:.   $z_1\leftarrow z$                              ▹ Save for following calculations 11:.   $z\leftarrow z+\ell ·\cos \theta$ 12:.   if $z\le 0$ then                        ▹ Check if photon hits lower boundary 13:.     $\ell \leftarrow -\frac{z_1}{\cos \theta }$ 14:.     $z\leftarrow 0$ 15:.     $reflected\leftarrow \mathbf{true}$ 16:.   end if 17:.   if $z\ge L$ then                        ▹ Check if photon hits upper boundary 18:.     $\ell \leftarrow \frac{L-z_1}{\cos \theta }$ 19:.     $z\leftarrow L$ 20:.     $reflected\leftarrow \mathbf{true}$ 21:.   end if 22:.   CalcRadContrib($z_1$, z, $\cos \theta$, $\sin \theta$, $\cos \varphi$, $\sin \varphi$, w, ℓ) 23:.   $w\leftarrow w·e^{-{\mu }_a\ell }$ 24:.   if reflected then                         ▹ Reflect photon and reduce weight 25:.     $\cos \theta \leftarrow -\cos \theta$ 26:.     $w\leftarrow w·CalcReflec\left(\theta \right)$ 27:.     $reflected\leftarrow \mathbf{false}$ 28:.   else                                        ▹ Scattering 29:.     ScatPhot($\cos \theta$, $\sin \theta$, $\cos \varphi$, $\sin \varphi$) 30:.   end if 31:.   if $w<ϵ$ then 32:.     RussianRoulette(q)           ▹ Break out of while loop if terminated by Russian roulette 33:.   end if 34:. end while

Algorithm 2 Calculate radiance contributions using LE

1:. function CalcRadContrib($z_1$, z, $\cos \theta$, $\sin \theta$, $\cos \varphi$, $\sin \varphi$, w, ℓ)   2:.   if $z\le 0$ or $z\ge L$ then                          ▹ Exclude direct detector hits   3:.     return   4:.   end if   5:.   for $i\leftarrow 1,n$ do                                 ▹ Loop over ${\theta }_i$ values   6:.     Rf ← CalcReflec$\left({\theta }_i\right)$   7:.     $e\leftarrow \exp \left(-{\mu }_{t}L/\left|\cos {\theta }_i\right|\right)$   8:.     for $j\leftarrow 1,m$ do                               ▹ Loop over ${\varphi }_j$ values   9:.       ${\widehat{s}}_p·\widehat{s}\leftarrow \sin \theta \sin {\theta }_i(\cos \varphi \cos {\varphi }_j+\sin \varphi \sin {\varphi }_j)+\cos \theta \cos {\theta }_i$ 10:.       ${\widehat{s}}_p·\tilde{s}\leftarrow \sin \theta \sin {\theta }_i(\cos \varphi \cos {\varphi }_j+\sin \varphi \sin {\varphi }_j)-\cos \theta \cos {\theta }_i$ 11:.       $w_s\leftarrow w·p\left({\widehat{s}}_p·\widehat{s}\right)·\exp \left(-{\mu }_t(L-z)/\left|\cos {\theta }_i\right|-{\mu }_a\ell \right)$          ▹ Equations (3) and (7) 12:.       ${\tilde{w}}_s\leftarrow w·p\left({\widehat{s}}_p·\tilde{s}\right)·\exp \left(-{\mu }_{t}z/\left|\cos {\theta }_i\right|-{\mu }_a\ell \right)$ 13:.       $I_{t_{i,j}}\leftarrow I_{t_{i,j}}+{\displaystyle \frac{1}{|\cos {\theta }_i|}}\left({\displaystyle \frac{w_s+{\tilde{w}}_s{R}_{f}e}{1-R_f^2{e}^2}}\right)$         ▹ Rad. contrib. on transm. side (Equation (8)) 14:.       $I_{r_{i,j}}\leftarrow I_{r_{i,j}}+{\displaystyle \frac{1}{|\cos {\theta }_i|}}\left({\displaystyle \frac{{\tilde{w}}_s+w_s{R}_{f}e}{1-R_f^2{e}^2}}\right)$         ▹ Rad. contrib. on reflect. side (Equation (9)) 15:.     end for 16:.   end for 17:. end function

Algorithm 3 Calculate radiance contributions using integral LE

1:. function CalcRadContrib($z_1$, z, $\cos \theta$, $\sin \theta$, $\cos \varphi$, $\sin \varphi$, w, ℓ)   2:.   for $i\leftarrow 1,n$ do                                   ▹ Loop over ${\theta }_i$ values   3:.     Rf ← CalcReflec$\left({\theta }_i\right)$   4:.     $e\leftarrow \exp \left(-{\mu }_t{\displaystyle \frac{L}{\cos {\theta }_i}}\right)$   5:.     for $j\leftarrow 1,m$ do                                 ▹ Loop over ${\varphi }_j$ values   6:.       ${\widehat{s}}_p·\widehat{s}\leftarrow \sin \theta \sin {\theta }_i(\cos \varphi \cos {\varphi }_j+\sin \varphi \sin {\varphi }_j)+\cos \theta \cos {\theta }_i$   7:.       ${\widehat{s}}_p·\tilde{s}\leftarrow \sin \theta \sin {\theta }_i(\cos \varphi \cos {\varphi }_j+\sin \varphi \sin {\varphi }_j)-\cos \theta \cos {\theta }_i$   8:.       $d_x\leftarrow {\mu }_a-{\mu }_t{\displaystyle \frac{\cos \theta }{\cos {\theta }_i}}$                    ▹ Denominators (Equations (29) and (30))   9:.       $d_y\leftarrow {\mu }_a+{\mu }_t{\displaystyle \frac{\cos \theta }{\cos {\theta }_i}}$ 10:.       if $d_x<{ϵ}_1$ then                       ▹ Check whether to use expansion for X 11:.         $X\leftarrow {\displaystyle \frac{1}{\cos {\theta }_i}}\exp \left({\mu }_t{\displaystyle \frac{z_1}{\cos {\theta }_i}}\right)\left(\ell -\frac{{\ell }^2{d}_x}{2}+\frac{{\ell }^3{d}_x^2}{6}\right)$                 ▹ Equation (31) 12:.       else 13:.         $X\leftarrow {\displaystyle \frac{1}{d_x\cos {\theta }_i}}\left(\exp \left\{{\mu }_t{\displaystyle \frac{z_1}{\cos {\theta }_i}}\right\}-\exp \left\{-{\mu }_a\ell +{\mu }_t{\displaystyle \frac{z}{\cos {\theta }_i}}\right\}\right)$    ▹ Equation (25) 14:.       end if 15:.       if $d_y<{ϵ}_2$ then                       ▹ Check whether to use expansion for Y 16:.         $Y\leftarrow {\displaystyle \frac{1}{\cos {\theta }_i}}\exp \left({\mu }_t{\displaystyle \frac{L-z_1}{\cos {\theta }_i}}\right)\left(\ell -\frac{{\ell }^2{d}_y}{2}+\frac{{\ell }^3{d}_y^2}{6}\right)$                 ▹ Equation (32) 17:.       else 18:.         $Y\leftarrow {\displaystyle \frac{1}{d_y\cos {\theta }_i}}\left(\exp \left\{{\mu }_t{\displaystyle \frac{L-z_1}{\cos {\theta }_i}}\right\}-\exp \left\{-{\mu }_a\ell +{\mu }_t{\displaystyle \frac{L-z}{\cos {\theta }_i}}\right\}\right)$    ▹ Equation (26) 19:.       end if 20:.       $r_{\mathrm{fac}}\leftarrow {\displaystyle \frac{R_f}{1-R_f^2{e}^2}}$ 21:.       $c_1\leftarrow p\left({\widehat{s}}_p·\widehat{s}\right)X+r_{\mathrm{fac}}\left(R_{f}p\left({\widehat{s}}_p·\widehat{s}\right)e^{2}X+p\left({\widehat{s}}_p·\tilde{s}\right)eY\right)$ 22:.       $c_2\leftarrow p\left({\widehat{s}}_p·\tilde{s}\right)Y+r_{\mathrm{fac}}\left(R_{f}p\left({\widehat{s}}_p·\tilde{s}\right)e^{2}Y+p\left({\widehat{s}}_p·\widehat{s}\right)eX\right)$ 23:.       $I_{r_{i,j}}\leftarrow I_{r_{i,j}}-{\mu }_{s}w{c}_1$                   ▹ Rad. contrib. on reflection side (Equation (27)) 24:.       $I_{t_{i,j}}\leftarrow I_{t_{i,j}}-{\mu }_{s}w{c}_2$                 ▹ Rad. contrib. on transmission side (Equation (28)) 25:.     end for 26:.   end for 27:. end function

3. Results and Discussion

The effectiveness of the iLE method was tested for different optical properties and compared with the LE method as a reference. To create representative comparison conditions, both methods have been integrated into the same self-developed MC framework, as described in Section 2. For prior validation of the MC code, the simulated results were first compared to an analytical solution of the RTE obtained using the $P_N$ method [24], and excellent agreement was found. In the context of the VRM benchmark test, the mean and variance of the simulated radiance were determined from 100 runs at a fixed computation time of 1.2 s per run, calculating 512 detection angles simultaneously. The Russian roulette parameters were set to $ϵ={10}^{-5}$ and $q=0.1$ for all test cases. All calculations were performed on an NVIDIA A10 graphics card using a thread count of 8192.

First, we consider the case of a thin layer with a thickness of $L=0.1\mathrm{m}\mathrm{m}$, a reduced scattering coefficient of ${\mu }_s^{\prime }={\mu }_s(1-g)=0.1{\mathrm{m}\mathrm{m}}^{-1}$, and an absorption coefficient of ${\mu }_a=0.1{\mathrm{m}\mathrm{m}}^{-1}$ (see Figure 2). The Henyey–Greenstein phase function [25] with an anisotropy factor of $g=0.8$ was chosen as the scattering function. In such a scenario, only a limited number of scattering events can be expected until the exit from the medium, so the iLE method, which, in contrast to the LE method, not only contributes to the radiance for individual scattering events but also allows the calculation of radiance contributions along the entire path segments, should be advantageous. This is confirmed by considering the variance in Figure 2. Above the total reflectance range $\left|\cos \theta \right|>\left|\cos {\theta }_c\right|$, the iLE method provides a variance improvement in reflectance of up to one order of magnitude and in transmittance of up to three orders of magnitude. It should be noted that the time dependence for both variance reduction methods was examined for all examples considered here and it was found, as expected, that the variance of both methods averaged over all angles decreases approximately with $1/N$, where N is the number of simulated photons. Consequently, the variance improvement corresponds to a reduction in simulation time of the same order of magnitude. As there tend to be very few interaction points up to the point of exit, especially in the transmission direction, this is also where the greatest improvement over the LE method is observed. Of particular interest is the range around $\left|\cos \theta \right|\approx 0.93$, which corresponds to the cosine of the angle of incidence ${\theta }_0^{\prime }$ after initial refraction at the interface. In this angular range, a peak in the mean radiance can be observed in both the refraction and transmission directions. The incident light angle of 30° combined with the phase function that favours forward scattering leads to photons of lower scattering order inducing an extended maximum, which can be mapped more effectively by the iLE method. It is important to note that the radiance is measured directly below the surface. Thus, for example, when comparing the simulated data with measured reflectance curves, only the region above the cosine of the critical angle of total reflection would be relevant since the radiance within the total reflection region would not be visible from the outside. Nevertheless, in the area of total reflection, the two methods show an equalisation, with the LE method being even slightly better for most of the angular range. However, for shallow viewing angles, which are almost parallel to the surface, the trend is reversed, as the long detection paths and the associated enormous attenuation mean that the individual LE weights make hardly any contribution, so the iLE method again proves to be superior by several orders of magnitude. In general, a further variance improvement via the iLE method compared to LE can be achieved by fine-tuning the parameters of the Russian roulette. In particular, the $ϵ$ threshold can significantly reduce the variance for computationally intensive iLE weights, as it prevents the propagation of photons that do not contribute significantly to the radiance.

For the next comparison, we examined the same parameters as those depicted in Figure 2, but this time with isotropic scattering, where $g=0$. This causes two pivotal changes within the mean radiance, as shown in Figure 3. First, the radiance curves in the transmission and reflection directions are now almost identical, since isotropic “sources” are almost evenly distributed across the thin layer. Secondly, the once visible peak around $\left|\cos \theta \right|\approx 0.93$ disappears. Moreover, reducing the g-factor also leads to a reduction in ${\mu }_s$, which further reduces the likelihood of scattering events, now giving the iLE method a lower variance across the full angular range. Note that since in this instance, the mean transmitted and reflected radiances are nearly identical, the variance of each method’s transmission and reflectance curve are also superimposed.

For a final test, see Figure 4, we again examine the same parameter set as in the first comparison, but for a substantially higher layer thickness of $L=10.0\mathrm{m}\mathrm{m}$ and a lower absorption of ${\mu }_a=0.01{\mathrm{m}\mathrm{m}}^{-1}$. Consequently, we anticipate a significant rise in potential interaction points in the medium and a marked decrease in attenuation along the propagation paths. This results in a noteworthy enhancement in the statistics for the LE technique, leading to almost similar variance for both VRMs across the majority of the angular range. Nonetheless, despite this specific situation, the iLE method still attains a substantial improvement of several orders of magnitude in variance for the relatively flat viewing angles.

4. Conclusions

A novel VRM, based on the integral form of the RTE, was introduced and benchmarked against the well-established LE method as a reference. To this end, the iLE and LE techniques were integrated into a self-developed Monte Carlo framework to conduct a thorough assessment of both methods’ statistical efficacy. The fundamental advantage of the iLE method lies in its integral nature, in which, unlike the LE technique, not only do individual scattering events contribute to the total radiance, but each individual path element adds integrally to the overall signal. This is also reflected in the test cases that have been examined. For a thin layer of 0.1 mm, where only a limited number of scattering events occur within the medium, a variance enhancement of up to three orders of magnitude could be observed below the critical angle of total reflection. This translates to a reduction in the simulation time by the same order of magnitude. Furthermore, for nearly horizontal viewing angles, where the attenuation along the detection path is exceptionally high, the variance could be improved even further for all scenarios under study. However, an increased layer thickness of 10 mm results in a substantial improvement in the statistics of the LE methods, so that in this case, both approaches yield comparable results, except for at shallow viewing angles, where the iLE method continues to provide significantly lower variance. In conclusion, the iLE method can be seen as an enhanced LE that allows a major improvement in statistics for media with low optical density and when near horizontal viewing angles are of interest, while still providing similarly good results for all other scenarios. Beyond the apparent application of simulating AR and AT data for arbitrary sets of optical parameters, the method is also predestined for solving inverse problems, since the simulated results yield perfectly smooth curves, owing to the fact that each path segment contributes to the full angular range. In forthcoming studies, the authors will strive to expand the methodology to structured illumination to enable its application in spatial frequency domain imaging (SFDI). In this domain, a considerable advantage of the introduced method is expected especially for large spatial frequencies, for which conventional Monte Carlo codes are very inefficient. Furthermore, there is a clear interest in extending the iLE method for spatially resolved radiance calculations to increase its general applicability.

Author Contributions

Conceptualization, D.H. and D.R.; methodology, D.H., D.R. and A.L.; software, D.H. and D.R.; validation, D.H. and D.R.; formal analysis, D.H. and D.R.; investigation, D.H. and D.R.; resources, A.K.; data curation, D.H. and D.R.; writing—original draft preparation, D.H.; writing—review and editing, D.R. and A.K.; visualization, D.H.; supervision, A.K.; project administration, D.H. and D.R.; funding acquisition, A.K. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All data sets are synthetic and produced by the published MC program. Archived source code at the time of publication is available as a Zenodo snapshot (https://doi.org/10.5281/zenodo.10300609, accessed on 20 December 2023).

Acknowledgments

The authors would like to thank Lars Heiler for reviewing the theory section and the corresponding GPU code.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Schematic representation showing the composition of local estimate weights for a single scattering event and detection direction indicating the radiance contributions at the respective interface. To calculate the transmittance or reflectance in the detection direction [Forumla omitted. See PDF.] after refraction, the last interaction at the interface must be taken into account according to Equation (10) or Equation (11), respectively. The photons are launched from a collimated source at an angle [Forumla omitted. See PDF.] onto the lower boundary surface ([Forumla omitted. See PDF.]).

Figure 2. Mean and variance of the simulated radiance in reflection (R) and transmission (T) for the parameters [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. For the externally measurable reflectance or transmittance, only the angular range to the right of the grey dotted line, which marks the critical angle of total reflection, is of relevance. The radiance shown here can be converted to the AT or AR using Equation (10) or Equation (11), respectively.

Figure 3. Mean and variance of the simulated radiance in reflection (R) and transmission (T) for the parameters [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. For the externally measurable reflectance or transmittance, only the angular range to the right of the grey dotted line, which marks the critical angle of total reflection, is of relevance. The radiance shown here can be converted to the AT or AR using Equation (10) or Equation (11), respectively.

Figure 4. Mean and variance of the simulated radiance in reflection (R) and transmission (T) for the parameters [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. For the externally measurable reflectance or transmittance, only the angular range to the right of the grey dotted line, which marks the critical angle of total reflection, is of relevance. The radiance shown here can be converted to the AT or AR using Equation (10) or Equation (11), respectively.

Appendix A. Derivation of the Integral Kernel

To derive the kernel function (19), we have to solve the radiative transfer Equation (14) under the reflecting boundary conditions (15) and (16) for a slab of thickness L. We introduce the parametrisation (A1)$$\widehat{K}(\ell ,\widehat{s}):=K(z+\ell \mu ,\widehat{s}).$$ In particular, we have $\widehat{K}(0,\widehat{s})=K(z,\widehat{s})$, yielding with the desired radiance. Inserting $\widehat{K}$ into (14) results in (A2)$$e^{-{\mu }_t\ell }{\partial }_{\ell }\left(\widehat{K}{e}^{{\mu }_t\ell }\right)=p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)\delta (z-z^{\prime }+\ell \mu ).$$ In the following, we consider two cases. First, for $\mu >0$, we obtain from $z+\ell \mu =0$ the value ${\ell }_1=-z/\mu <0$. Integrating (A2) over $({\ell }_1,0)$ gives (A3)$$K(z,\widehat{s})=K(0,\widehat{s})\exp \left(-{\mu }_t\frac{z}{\mu }\right)+p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)G(z,\mu ),\mu >0,$$ where (A4)$$G(z,\mu ):=\frac{1}{\left|\mu \right|}\exp \left(-{\mu }_t\frac{z-z^{\prime }}{\mu }\right)\Theta \left(\frac{z-z^{\prime }}{\mu }\right).$$ Second, for $\mu <0$, setting $z+\ell \mu =L$ gives ${\ell }_2=(L-z)/\mu <0$. Again, integrating (A2) over $({\ell }_2,0)$ leads to (A5)$$K(z,\widehat{s})=K(L,\widehat{s})\exp \left({\mu }_t\frac{L-z}{\mu }\right)+p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)G(z,\mu ),\mu <0.$$ We note that the values $K(0,\widehat{s})$ in (A3) and $K(L,\widehat{s})$ in (A5) are not yet known. To incorporate the boundary conditions (15) and (16), we additionally need the values for the reflected radiance $K(0,-\mu ,\varphi )$ and $K(L,\mu ,\varphi )$. From (A5), we find (A6)$$K(0,\mu ,\varphi )=K(L,\mu ,\varphi )\exp \left({\mu }_t\frac{L}{\mu }\right)+\frac{p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)}{\left|\mu \right|}\exp \left({\mu }_t\frac{z^{\prime }}{\mu }\right),\mu <0,$$ which is the same as (A7)$$K(0,-\mu ,\varphi )=K(L,-\mu ,\varphi )\exp \left(-{\mu }_t\frac{L}{\mu }\right)+\frac{p\left(\tilde{s}·{\widehat{s}}^{\prime }\right)}{\left|\mu \right|}\exp \left(-{\mu }_t\frac{z^{\prime }}{\mu }\right),\mu >0,$$ with $\tilde{s}:={(\sqrt{1-{\mu }^2}\cos \varphi ,\sqrt{1-{\mu }^2}\sin \varphi ,-\mu )}^T$. Similarly, we obtain from (A3) the boundary value (A8)$$K(L,\mu ,\varphi )=K(0,\mu ,\varphi )\exp \left(-{\mu }_t\frac{L}{\mu }\right)+\frac{p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)}{\mu }\exp \left(-{\mu }_t\frac{L-z^{\prime }}{\mu }\right),\mu >0.$$ We now can use the BC to obtain the following system of linear equations (A9)$$\begin{array}{cc}\hfill & K(0,-\mu ,\varphi )=R_f\left(\mu \right)K(L,\mu ,\varphi )\exp \left(-{\mu }_t\frac{L}{\mu }\right)+\frac{p\left(\tilde{s}·{\widehat{s}}^{\prime }\right)}{\mu }\exp \left(-{\mu }_t\frac{z^{\prime }}{\mu }\right),\hfill \end{array}$$ (A10)$$\begin{array}{cc}\hfill & K(L,\mu ,\varphi )=R_f\left(\mu \right)K(0,-\mu ,\varphi )\exp \left(-{\mu }_t\frac{L}{\mu }\right)+\frac{p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)}{\mu }\exp \left(-{\mu }_t\frac{L-z^{\prime }}{\mu }\right),\hfill \end{array}$$ where $\mu >0$ in both equations, for the unknowns $K(0,-\mu ,\varphi )$ and $K(L,\mu ,\varphi )$. The corresponding solution is given by (A11)$$\begin{array}{cc}\hfill & K(0,-\mu ,\varphi )=\frac{1}{\mu }\frac{p\left(\tilde{s}·{\widehat{s}}^{\prime }\right)e^{-{\mu }_t{z}^{\prime }/\mu }+R_f\left(\mu \right)p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)e^{{\mu }_t(z^{\prime }-2L)/\mu }}{1-R_f^2\left(\mu \right)e^{-2{\mu }_{t}L/\mu }},\hfill \end{array}$$ (A12)$$\begin{array}{cc}\hfill & K(L,\mu ,\varphi )=\frac{1}{\mu }\frac{p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)e^{{\mu }_t(z^{\prime }-L)/\mu }+R_f\left(\mu \right)p\left(\tilde{s}·{\widehat{s}}^{\prime }\right)e^{-{\mu }_t(z^{\prime }+L)/\mu }}{1-R_f^2\left(\mu \right)e^{-2{\mu }_{t}L/\mu }}.\hfill \end{array}$$ With these values at hand, we can compute $K(0,\mu ,\varphi )=R_f\left(\mu \right)K(0,-\mu ,\varphi )$ ($\mu >0$) as well as $K(L,\mu ,\varphi )=R_f(-\mu )K(L,-\mu ,\varphi )$ ($\mu <0$), yielding in view of (A3) and (A5), (A13)$$\begin{array}{cc}\hfill & K(z,\widehat{s})=R_f\left(\mu \right)K(0,-\mu ,\varphi )\exp \left(-{\mu }_t\frac{z}{\mu }\right)+p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)G(z,\mu ),\mu >0,\hfill \end{array}$$ (A14)$$\begin{array}{cc}\hfill & K(z,\widehat{s})=R_f(-\mu )K(L,-\mu ,\varphi )\exp \left({\mu }_t\frac{L-z}{\mu }\right)+p\left(\widehat{s}·{\widehat{s}}^{\prime }\right)G(z,\mu ),\mu <0.\hfill \end{array}$$ Then, inserting (A4) and combining the expressions for the two cases results in the final form, see Equation (19), for the kernel function.

References

1. Vincendon, M.; Langevin, Y.; Poulet, F.; Bibring, J.P.; Gondet, B. Recovery of surface reflectance spectra and evaluation of the optical depth of aerosols in the near-IR using a Monte Carlo approach: Application to the OMEGA observations of high-latitude regions of Mars. J. Geophys. Res. Planets; 2007; 112, E08S13. [DOI: https://dx.doi.org/10.1029/2006JE002845]

2. Barker, H.; Goldstein, R.; Stevens, D. Monte Carlo Simulation of Solar Reflectances for Cloudy Atmospheres. J. Atmos. Sci.; 2003; 60, pp. 1881-1894. [DOI: https://dx.doi.org/10.1175/1520-0469(2003)060<1881:MCSOSR>2.0.CO;2]

3. Bréon, F.M. Reflectance of Broken Cloud Fields: Simulation and Parameterization. J. Atmos. Sci.; 1992; 49, pp. 1221-1232. [DOI: https://dx.doi.org/10.1175/1520-0469(1992)049<1221:ROBCFS>2.0.CO;2]

4. McKee, T.; Cox, S. Simulated Radiance Patterns for Finite Cubic Clouds. J. Atmos. Sci.; 1976; 33, pp. 2014-2020. [DOI: https://dx.doi.org/10.1175/1520-0469(1976)033<2014:SRPFFC>2.0.CO;2]

5. Janecek, M.; Moses, W. Simulating Scintillator Light Collection Using Measured Optical Reflectance. IEEE Trans. Nucl. Sci.; 2010; 57, pp. 964-970. [DOI: https://dx.doi.org/10.1109/TNS.2010.2042731]

6. Janecek, M.; Moses, W. Optical Reflectance Measurements for Commonly Used Reflectors. IEEE Trans. Nucl. Sci.; 2008; 55, pp. 2432-2437. [DOI: https://dx.doi.org/10.1109/TNS.2008.2001408]

7. Trigila, C.; Moghe, E.; Roncali, E. Technical Note: Standalone application to generate custom reflectance Look-Up Table for advanced optical Monte Carlo simulation in GATE/Geant4. Med. Phys.; 2021; 48, pp. 2800-2808. [DOI: https://dx.doi.org/10.1002/mp.14863] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33772816]

8. Barajas, O.; Ballangrud, A.; Miller, G.; Moore, R.; Tulip, J. Monte Carlo modelling of angular radiance in tissue phantoms and human prostate: PDT light dosimetry. Phys. Med. Biol.; 1997; 42, 1675. [DOI: https://dx.doi.org/10.1088/0031-9155/42/9/001]

9. Disney, M.; Lewis, P.; North, P. Monte Carlo ray tracing in optical canopy reflectance modelling. Remote Sens. Rev.; 2000; 18, pp. 163-196. [DOI: https://dx.doi.org/10.1080/02757250009532389]

10. North, P. Three-dimensional forest light interaction model using a Monte Carlo method. IEEE Trans. Geosci. Remote Sens.; 1996; 34, pp. 946-956. [DOI: https://dx.doi.org/10.1109/36.508411]

11. Cooper, K.; Smith, J. A Monte Carlo Reflectance Model for Soil Surfaces with Three-Dimensional Structure. IEEE Trans. Geosci. Remote Sens.; 1985; GE-23, pp. 668-673. [DOI: https://dx.doi.org/10.1109/TGRS.1985.289385]

12. Kurt, M.; Edwards, D. A Survey of BRDF Models for Computer Graphics. SIGGRAPH Comput. Graph.; 2009; 43, pp. 1-7. [DOI: https://dx.doi.org/10.1145/1629216.1629222]

13. Guarnera, D.; Guarnera, G.; Ghosh, A.; Denk, C.; Glencross, M. BRDF Representation and Acquisition. Comput. Graph. Forum; 2016; 35, pp. 625-650. [DOI: https://dx.doi.org/10.1111/cgf.12867]

14. Flock, S.; Patterson, M.; Wilson, B.; Wyman, D. Monte Carlo modeling of light propagation in highly scattering tissues. I. Model predictions and comparison with diffusion theory. IEEE Trans. Biomed. Eng.; 1989; 36, pp. 1162-1168. [DOI: https://dx.doi.org/10.1109/TBME.1989.1173624]

15. Wong, B.; Mengüç, M. Comparison of Monte Carlo techniques to predict the propagation of a collimated beam in participating media. Numer. Heat Transf. Part B Fundam.; 2002; 42, pp. 119-140. [DOI: https://dx.doi.org/10.1080/10407790190053860]

16. Wang, L.; Jacques, S.; Zheng, L. MCML—Monte Carlo modeling of light transport in multi-layered tissues. Comput. Methods Progr. Biomed.; 1995; 47, pp. 131-146. [DOI: https://dx.doi.org/10.1016/0169-2607(95)01640-F]

17. Pincus, R.; Evans, K. Computational Cost and Accuracy in Calculating Three-Dimensional Radiative Transfer: Results for New Implementations of Monte Carlo and SHDOM. J. Atmos. Sci.; 2009; 66, pp. 3131-3146. [DOI: https://dx.doi.org/10.1175/2009JAS3137.1]

18. Iwabuchi, H. Efficient Monte Carlo Methods for Radiative Transfer Modeling. J. Atmos. Sci.; 2006; 63, pp. 2324-2339. [DOI: https://dx.doi.org/10.1175/JAS3755.1]

19. Buras, R.; Mayer, B. Efficient unbiased variance reduction techniques for Monte Carlo simulations of radiative transfer in cloudy atmospheres: The solution. J. Quant. Spectrosc. Radiat. Transf.; 2011; 112, pp. 434-447. [DOI: https://dx.doi.org/10.1016/j.jqsrt.2010.10.005]

20. Marchuk, G.; Mikhailov, G.; Nazareliev, M.; Darbinjan, R.; Kargin, B.; Elepov, B. The Monte Carlo Methods in Atmospheric Optics; Springer: Berlin/Heidelberg, Germany, 1980; 208.

21. Wilson, B.; Adam, G. A Monte Carlo model for the absorption and flux distributions of light in tissue. Med. Phys.; 1983; 10, pp. 824-830. [DOI: https://dx.doi.org/10.1118/1.595361]

22. Martelli, F.; Binzoni, T.; Bianco, S.D.; Liemert, A.; Kienle, A. Light Propagation through Biological Tissue and Other Diffusive Media: Theory, Solutions, and Validations; 2nd ed. SPIE: Bellingham, WA, USA, 2022; 51.

23. Machida, M. Rotated Reference Frames in Radiative Transport Theory; Progress in Optics Elsevier: Amsterdam, The Netherlands, 2023; [DOI: https://dx.doi.org/10.1016/bs.po.2023.06.001]

24. Liemert, A.; Reitzle, D.; Kienle, A. Hybrid method for solving the radiative transport equation. Proceedings of the Diffuse Optical Spectroscopy and Imaging IX; Munich, Germany, 25–29 June 2023; Contini, D.; Hoshi, Y.; O’Sullivan, T.D. International Society for Optics and Photonics, SPIE: Bellingham, WA, USA, 2023; Volume 12628, 126281G. [DOI: https://dx.doi.org/10.1117/12.2670672]

25. Henyey, L.; Greenstein, J. Diffuse radiation in the galaxy. Astrophys. J.; 1941; 93, pp. 70-83. [DOI: https://dx.doi.org/10.1086/144246]