Incorporating reflection boundary conditions in the Neumann series radiative transport equation: application to photon propagation and reconstruction in diffuse optical imaging

: We propose a formalism to incorporate boundary conditions in a Neumann-series-based radiative transport equation. The formalism accurately models the reﬂection of photons at the tissue-external medium interface using Fresnel’s equations. The formalism was used to develop a gradient descent-based image reconstruction technique. The proposed methods were implemented for 3D diﬀuse optical imaging. In computational studies, it was observed that the average root-mean-square error (RMSE) for the output images and the estimated absorption coeﬃcients reduced by 38% and 84%, respectively, when the reﬂection boundary conditions were incorporated. These results demonstrate the importance of incorporating boundary conditions that model the reﬂection of photons at the tissue-external medium interface.


Introduction
In optical imaging modalities such as diffuse optical imaging (DOI) and fluorescence imaging, using the detector measurements to estimate the optical coefficients of the imaged tissue requires an accurate model for photon propagation. In these modalities, the detectors are usually placed on the surface of the tissue leading to a refractive index mismatch as photons exit from the tissue to the external medium. Due to this mismatch, a certain percentage of photons are reflected back into the tissue. Accounting for the boundary conditions arising due to this refractive index mismatch is necessary for accurate modeling of photon propagation [1,2]. It has been observed that not accounting for boundary conditions accurately leads to erroneous modeling of scattered light [3], and 50% or more errors in estimating the optical coefficients of the tissue [1].
There has been considerable research on incorporating boundary conditions while modeling photon transport in biological tissue, especially when using the radiative transport equation (RTE). Here we summarize this research with a focus on the boundary conditions implemented in existing methods. Analytical solutions for RTE have been studied where a Marshak-type boundary condition has been proposed. This condition incorporates reflection of photons accurately using Fresnel's equations [4][5][6]. However, analytical methods typically work only for specific geometries and thus numerical solutions to RTE have been widely explored. Using the complete form of the RTE to model photon transport numerically is computationally expensive [7]. Thus, typically, an approximate form of the RTE is used, which assumes that photons propagate diffusively through the tissue. These diffusion approximation (DA)-based methods are computationally faster, but have several limitations. For example, in media that have low scattering, high absorption or small geometry, the DA-based methods are inaccurate [8][9][10][11]. Consequently, these methods have limitations in describing photon transport in highly absorbing regions such as haematomas, void-like spaces such as ventricles and the subarachnoid-space, and for small-animal imaging [10][11][12][13][14]. Further, modeling the reflection of photons is non-trivial using the DA-based methods because while reflection is highly directional, the DA-based methods assume that photons propagate diffusively. While different boundary conditions including partial-current, extrapolated, and vacuum boundary conditions have been proposed, these have varying degrees of inaccuracy [1,15,16]. To correct for these inaccuracies, a corrected diffusion approximation method has been proposed [17,18]. However, given the limitations of DA-based methods, there is an important need for more accurate computational photon-transport models.
To model photon transport accurately, numerical methods that use the differential form of the RTE have been widely investigated. These methods represent important advances and have yielded promising results. However, these methods require solving many coupled equations leading to large computational times. For example, higher-order spherical harmonics-based methods (P N ) [19][20][21][22]  More recently boundary conditions that model the reflection of photons were implemented for the S N methods [23,30]. To overcome the computational complexity of the S N and P N methods, a simplified spherical harmonics (SP N ) method has been proposed [10,14,31,32]. The SP N methods implement partial-reflection boundary conditions. However, the SP N solutions also have several limitations [10].
To overcome the issue of solving multiple coupled equations, integral forms of the RTE have been explored. More specifically, a Neumann-series-based RTE was developed to model photon transport for both homogeneous [33] and heterogeneous scattering media [34]. This method requires solving only a single equation. Further, the method also provided improved accuracy for imaging setups where the DA-based methods have limitations [33]. Additionally, the method provided a novel intuitive framework to physically interpret the RTE solution. Each term in the Neumann-series RTE describes contributions from photons that have scattered a specific number of times. Thus, the method can numerically quantify the contribution from photons that, for example, do not scatter (ballistic photons), or scatter only a certain number of times. However, the existing Neumann-series method inaccurately assumes that there is no reflection of photons at the boundary, which, as mentioned above, can lead to inaccurate modeling of photon transport and inaccurate estimation of optical coefficients. Further, inaccurate modeling of boundary conditions also limits the experimental validation of the Neumann-series method.
To overcome the above-mentioned issue, we propose a novel Neumann-series RTE that models the reflection of photons at the boundary between tissue and external medium. We use this Neumann series RTE to investigate the improvement in estimating optical coefficients when reflection of photons at the boundary is modeled. For this purpose, we propose a novel reconstruction technique that incorporates the Neumann-series formalism to estimate the optical coefficients of the tissue. All these methods are developed in the context of a three-dimensional (3D) DOI imaging system. Preliminary versions of this work have been presented previously [35-37].

The Neumann-series RTE modeling the boundary conditions
The fundamental radiometric quantity that we describe using the RTE is the photon distribution function, a quantity that is proportional to radiance and quantifies the density of photons at a particular location and in a particular direction. Denote the photon distribution function at location r = (x, y, z) in directionŝ by w(r,ŝ). In terms of photons, the photon distribution function w(r,ŝ)∆V∆Ω can be interpreted as the number of photons contained in volume ∆V centered on the 3D vector r and traveling in solid angle ∆Ω about directionŝ at a given time t. Denote the absorption and scattering coefficients at location r by µ a (r) and µ s (r), respectively, and define µ tot (r) = µ a (r) + µ s (r), where µ tot is also referred to as the transport coefficient. Also, denote the speed of light in the medium by c m . Let Ξ(r,ŝ) denote a monochromatic mono-energetic emission source, and let f (ŝ,ŝ ; r) denote the scattering phase function. For this source, the RTE can be written as [33] s · ∇w(r,ŝ) + µ tot w(r,ŝ) = 1 c m Ξ(r,ŝ) + µ s (r) ∫ 4π dΩ f (ŝ,ŝ ; r)w(r,ŝ ) .
One choice for scattering phase function in biological tissue is the Henyey-Greenstein function.
f (ŝ,ŝ ; r) = 1 4π where g ∈ [−1, 1] is anisotropy factor. The RTE can also be represented in an integral form in terms of a Neumann series. This series consists of scattering and attenuation operators, denoted Similarly where λ corresponds to the length variable. A schematic demonstrating the scattering and attenuation operators for the 3D DOI system considered in this manuscript is given in Fig. 1. Using these notations, the integral form of the RTE is given by [33] w = XΞ + XKw, where w is short-hand notation for w(r,ŝ). This equation can be alternatively written in a Neumann series form as follows: The equation has a very intuitive physical interpretation; the contribution of photons that have scattered n times is given by the term X(KX) n Ξ in this series.
To incorporate the boundary conditions into the Neumann series, we start with a first principles treatment of light propagation in tissue. Similar to Schweigher et al.
[38], we define the boundary of the tissue to be an infinitesimally thin layer where only the reflection event occurs. This thin boundary layer is illustrated in Fig. 2. Due to the refractive index mismatch, photons incident on the boundary are reflected within this thin layer. Therefore, the thin layer acts as a source of photon emission. Let R denote the reflection operator. An alternative way to think about the reflection operation is to consider it as a scattering operation, except that the scattering phase function is given by the laws of reflection. Either of these two interpretations, when modeled in Eq. 5, leads to the following form for the RTE: The operator R is defined using Fresnel's laws of reflection. Letŝ ,ŝ andn denote the angle of incidence, reflection, and normal to the boundary surface. In accordance with Fresnel's laws of reflection, where δ(. . .) is the Kronecker delta function and where b(r) = 0 is the equation of the boundary. Thus, δ(b(r)) denotes that the reflection operator acts only on photons in the infinitesimally thin boundary layer. For example, if the boundary is a plane normal ton and at a perpendicular distance p from the origin, then b(r) = r ·n − p. Further, R(ŝ) is reflectivity coefficient. Let n i and n t denote refractive index of incident medium and transmitted medium, respectively. Also let θ i , θ t , and θ c denote incident, transmitted, and critical angle, respectively. Then, Taking the terms involving the distribution function on the left side of Eq. 7 yields A formal solution of the above equation is given by the Neumann series: where I is the identity operator. Similar to the originally proposed Neumann-series RTE, the various terms in the above expansion have physical interpretations, as shown in Fig. 3. For example, the term XRXΞ represent the photons that are reflected at the boundary and subsequently transmitted back into the medium. The photons that reflect back into the medium and then get scattered are represented by the term XKXRXΞ. The term XRXKXΞ denotes photons that are scattered in the medium and subsequently reflected from a boundary surface.

Implementation
The Neumann-series formalism was implemented for a 3-D DOI setup as shown in Fig. 1. The output image was measured on a pixellated contact detector. For computational reasons, the effect of reflection was accounted using only three terms that contained the reflection operator, as follows: The implementation of the Neumann series without reflection boundary conditions was similar to a previous implementation [33]. The mathematical expressions for the three additional terms containing the reflection operator is derived below. The physical interpretation of these terms is shown in Fig. 3. The first term XRXΞ describes photons that are emitted from the source, pass through the medium, incident on the exit boundary, and then reflected back into the medium. In our DOI setup, the laser source emits a unidirectional beam along the optical axis, considered as the z axis. Denote the plane where the laser source is incident by z = 0 and the length of the medium along the z axis by L. Denote the profile of the beam along the (x, y) coordinates by h(x, y) and let α denote the intensity of the laser source. Then the source term is [33] For this source term, the expression for the XRXΞ term is derived in Appendix A and is The second term XKXRXΞ describes photons that scatter once after reflection before reaching the detector surface. Representing this term using the local basis functions in angle requires considerable memory, as exemplified previously [33]. To overcome this issue, this term is implemented in the spherical harmonic basis. The distribution function can be written in the spherical harmonic basis as follows: where Y lm (ŝ) denote the spherical harmonic basis functions and W lm (r) denotes the spherical harmonic coefficients, which are Denote the term XRXΞ, given by Eq. 14, by a new source term Ξ . Applying the scattering operator (Eq. 3) to this source term gives Denote the scattering operator in spherical harmonics by D and the source term Ξ by χ . The term [KΞ ](r,ŝ) is denoted in the spherical harmonic basis as [D χ ] lm (r). Using Eq. 16 yields The attenuation operator is then applied numerically on this term. The expression for the attenuation operator in the spherical harmonic basis has been derived previously [33,39]. The third term XRXKXΞ describes photons that scatter once before being reflected at boundary. This boundary is the set of all planes surrounding the medium, as shown in Fig. 2, and is thus defined as i δ(r ·n i − p i ), where the index i denotes the different planes,n i denotes the unit normal vector to the i th plane and p i denotes the distance of the i th plane from the origin along the normal vector. With this notation, the expression for the XRXKXΞ term, as derived in Appendix B, is where From the computed distribution function, the transmitted flux on the detector face for the DOI system is computed as follows. Denote the x and y dimensions of each pixel by ∆x and ∆y, respectively. Denote the transmitted flux detected by the m th pixel with center at r m by Φ m . Assuming that the distribution function over a pixel is approximately constant, we get [33] where T(ŝ) = 1 − R(ŝ) is transmission coefficient. As an example, the flux due to the XRXKXΞ term is derived in Appendix B.

Designing the Neumann-series-based reconstruction algorithm accounting for boundary conditions
Denote the absorption/scattering coefficient at location r by the function µ(r). The function is discretized using the spatial basis function φ n (r), as below: where µ n denote the optical coefficient (scattering or absorption) for the n th basis function. The spatial basis functions could be the commonly used voxels, but the method is generalizable to other spatial basis functions too. Denote the N-dimensional vector of the coefficients µ n by µ µ µ. Denote the image acquired by the pixelated detector in the DOI setup by the M-dimensional vector g. Our objective is to estimate µ µ µ given g. For this purpose, we derive a gradient-descent-based approach.
Denote the mean noiseless image as a function of the scattering and absorption coefficient bȳ g(µ µ µ). In our reconstruction approach, the objective is to estimate the value of µ µ µ that minimizes the L 2 norm of the error between the measured data g and theḡ(µ µ µ). Mathematically, the objective is to estimateμ µ µ that satisfiesμ where Ψ(µ µ µ) = g −ḡ(µ µ µ) 2 2 is the objective function, and where x 2 2 denotes the square of the L 2 norm of the vector x. To implement the gradient-descent method, we need to calculate the derivative of Ψ(µ µ µ), which can be written as The calculation of the above equation requires computing the derivative ofḡ(µ µ µ) with respect to µ a and µ s . Let h m (r,ŝ) denote the detector response function of the m th detector pixel. Then, the m th component ofḡ(µ µ µ), denoted byḡ m , is where ( , ) denotes the inner product of two vectors. Taking the derivative on both sides with respect to µ n yields From Eq. 7, As derived previously [40], terms containing partial derivatives of X and K are and where = 1 for µ = µ s and = 0 for µ = µ a , and the operator K 1 is defined as Further, as derived in Appendix C Substituting the expressions from Eqs. (29)-(33) in Eq. (28) and rearranging the terms yields From Eq. (7), XΞ + XRw + XKw = w. Thus Define S = −c m φ n w + φ n K 1 w. The above equation can then be rewritten as Comparison to Eq. 7 reveals that the expression for the gradient is the same as the original RTE, but with a different source term. Therefore, the gradient of the distribution function is computed by simply executing the Neumann-series RTE but with the source term S. From the gradient of the distribution function and using Eq. (27), the gradient of g m with respect to µ n is obtained. Using a simple iterative gradient-descent approach, the computed gradient is used to update the values of the optical coefficients in each iteration until convergence is achieved. Thus the absorption and scattering coefficients are estimated. Specifically, for a homogeneous medium, the N-dimensional coefficients vector µ µ µ in the above derivation becomes a scalar value µ. The calculation of ∂ḡ m ∂µ remains the same as the above derivation.

Experimental setup
The proposed Neumann-series RTE approach was validated for a 3D DOI setup (Fig. 4). The scattering medium was a cube with each side of length 2 cm. The scattering and absorption coefficients in the medium were 1 cm −1 and 0.01 cm −1 . The medium was small geometry, had a relatively low scattering coefficient, and a collimated source. The choice for this setup was made to study the performance of the proposed method in a scenario where the diffusion approximation-based methods are known to have limitations. In the simulation setup, across different experiments, the refractive index of the media was varied from 1.1 to 1.5, while the refractive index of the external medium was kept as 1. The anisotropy factor g for the scattering medium was set to 0. A pixellated 2D contact detector with 20 × 20 pixels acquired the output image.
Our first objective was to investigate the performance of the proposed Neumann-series method in modeling photon transport. For this purpose, output images and fluence fields were generated using the proposed method. These were compared to the images and the fields obtained with the . For the MC study, the medium was discretized into 40 × 40 × 40 voxels. As in previous studies, the output using the MC technique was considered as the gold standard [33,34]. The difference between the output obtained using the MC and the proposed Neumann-series technique over the entire image was quantified using the normalized root mean square error (RMSE) metric. Let M denote the number of pixels in the acquired image. In our setup, M = 400. Let g m, RT E and g m, MC denote the measurements acquired by the m th detector pixel. Then the normalized RMSE was defined as follows: We also examined the effect on accuracy of modeling photon transport when the reflection of photons at the boundary was not accounted. For this purpose, the output images were obtained using an existing Neumann-series method that assumed that the photons are completely absorbed at the boundary [33]. The RMSE for these output images were compared to those obtained using the proposed Neumann-series technique, thus quantifying the improvement in accuracy. The second objective was to examine the effect on the accuracy of the estimated optical coefficients when the reflection of photons at the boundary was modeled. For this purpose, first the image data for the DOI setup in Fig. 4 was generated using the MC technique. Next, the absorption coefficient of the medium was estimated from this generated image using two reconstruction methods; the proposed Neumann-series-based reconstruction method and a modified Neumannseries-based reconstruction method that did not contain the reflection operator [33]. The scattering coefficient was assumed to be known. For this setup, the spatial basis function (Eq. (23)) was simply the entire support of the object so that N = 1 and the objective was to only estimate a single value µ a . The RMSE of the estimated absorption coefficient using the two reconstruction methods were compared.

Results
Results from the experiments investigating the accuracy of photon transport are presented in Figs. 5 and 6. In Figs. 5a-d, the linear profile along the center of the output image obtained using the proposed Neumann series approach and using the MC technique are plotted for different amounts of refractive index mismatch. It is observed that the linear profiles overlap, demonstrating the accuracy of the proposed technique. The RMSE between the output images obtained using the MC and proposed Neumann-series RTE method is shown in Table 1. The corresponding RMSE using the existing Neumann-series technique that does not model the reflection of photons [33] is also shown. It is observed that the RMSE values were lower using the proposed Neumann series technique with the average RMSE reduced by 38%.
The fluence on a 2D plane defined at y = 0.9 cm using the MC and proposed Neumann-series RTE method is shown in Fig. 6. In the third column of Fig. 6, the fluence contours for the Neumann-series RTE and the MC method are overlaid on top of each other, facilitating a visual comparison. It is observed that for different values of refractive index mismatch, the contour fields using the MC and RTE methods lie approximately on top of each other.
Results from the experiments investigating the reconstruction method are presented in Fig. 7. It is observed that using the proposed reconstruction algorithm, the estimated absorption coefficient was close to the true absorption coefficient of 0.01 cm −1 for all values of refractive index, with an average normalized RMSE of 25.1%. However, on using the modified reconstruction algorithm that did not contain the reflection operator, the estimated absorption coefficient had a large error for all values of refractive index, with an average RMSE of 156.5%. On an average, we found that the mean RMSE of the estimated absorption coefficient reduced by 84%. Fig. 6. Comparison of the fluence calculated using the MC and proposed Neumann-series RTE-based methods. The refractive index of the medium is n=1.2 (first row), n=1.3 (second row), n=1.4 (third row), and n=1.5 (fourth row). In the first column (a, d, g, j), the logarithm of the amplitude of the fluence obtained with the MC method is plotted. In the second column (b, e, h, k), the logarithm of the amplitude of the fluence obtained with Neumann-series RTE is plotted. In the third column (c, f, i, l), the contour of the fluence obtained with the MC and Neumann-series RTE is plotted at 3 dB spacing.

Discussions
The first goal of this paper was to propose a Neumann-series RTE that modeled the reflection of photons at the boundary surface. The results in Fig. 5 and 6 demonstrate that the proposed Neumann-series method provided an accurate modeling of the photon transport for the considered experimental setup. Further, results in Table 1 show that the proposed technique yields improved accuracy in modeling photon transport in comparison to when the reflection of photons at the boundary is not modeled. Thus, these results demonstrate that accounting for the reflection of photons at the boundary is essential to model photon transport accurately. The second goal of this paper was to investigate the effect of modeling reflection at the boundary on the accuracy of the estimated optical coefficients. The results in Fig. 7 show that the optical coefficients estimated when the reflection was modeled were substantially more accurate than when the reflection at the boundary was not modeled. These results demonstrate the importance of modeling the reflection of photons at the boundary of the tissue and the external medium.
To implement the proposed Neumann-series method, we considered only three additional terms in comparison to the original Neumann series. In each of these terms, the reflection operator was present once, under the assumption that terms that contain two or more instances of the reflection operator would not have a significant contribution to the output. This is because if the reflection coefficient R is not very high such that R 2 << 1, terms that contain the reflection operator twice could be ignored. For example, for the experimental setup in this manuscript, the plot of the reflection coefficient as a function of the incident angle is plotted in Fig. 8. We observe that the value of R is close to 0 for a majority of the angles of incidence. Another reason for choosing these three terms was that computing the other terms containing the reflection operator would require implementing the reflection operator in spherical harmonic basis. This is complicated because the reflection operator is highly directional. The results in Fig. 5 show that the output using the MC and RTE approaches match well when only these three additional terms are used. Because the refractive index mismatch in the simulation study setup are representative of actual imaging setups, these results provide strong evidence for using only three additional terms in the Neumann-series RTE.
The Neumann-series method can handle different phase functions in a way similar to the other more conventional integro-differential methods to solve the RTE. This is because the scattering operator is in an integral form in both the Neumann series and the integro-differential form of the RTE. Also, the Neumann-series method method can handle complex illumination patterns. Note that the source term is modeled by the Ξ(r,ŝ) term in the Neumann-series RTE. The DOI setup used in this paper consisted of a collimated monochromatic laser source, so that the source term was defined as in Eq. (13). For a different illumination pattern, this term would need to be modified accordingly. Finally, while the setup considered in this manuscript assumed a time-independent source, the Neumann-series RTE is applicable to time-dependent sources. For this purpose, note that the RTE can be represented in the frequency domain for each frequency component ν by replacing µ tot in Eq. (1) by µ tot + iν cm , where i denotes the imaginary unit. Thus, the Neumann-series RTE can be solved for each frequency component by simply replacing µ tot by µ tot + iν cm [33]. Consequently, the proposed Neumann-series method could be used to model photon propagation in a wide range of imaging setups. We do note that for the different imaging setups, the number of terms required in the Neumann-series might be different and not known a priori. For this purpose, a criterion has been developed that enables assessing whether the Neumann-series RTE output has converged [33]. Using this criterion allows an adaptive determination of whether the Neumann series has converged or whether more terms need to be computed to obtain an accurate output.
For the experimental setup considered in this manuscript, the DA-based methods are inaccurate [8][9][10][11]33]. Our results demonstrate the accuracy and resultant advantage of the proposed Neumann series method in modeling photon transport in this setup. However, when the scattering coefficients are high, or the medium has a large geometry, the Neumann series method has large memory and computational requirements. Previously we have observed that the method is accurate and practical when µ s L is around 5 [33], where L denotes the length of the medium. This is a limitation of the proposed method. In these cases, the DA-based methods yield accurate results. Thus, the proposed method could be integrated with the DA methods to yield hybrid approaches, similar to the RTE-DA method [42] and MC-DA method [43] proposed previously. Further, rapid advances in the high-performance computing hardware technology are being observed. For example, a NVIDIA graphics processing units (GPU)-based implementation of the RTE was proposed with the NVIDIA C2050 GPUs in 2012 [34]. The current generation of GPUs, namely the NVIDIA Volta, have about fifteen-times the processing power and five-times the memory of the C2050. Since the challenges with the Neumann-series method are mainly computational, we anticipate that these advances in high-performance computing will alleviate these challenges in the near future.
The implementation of the Neumann series approach with boundary conditions enables an experimental validation of this technique. This experimental validation is an important direction of future research. Another important frontier is the implementation of this method for heterogeneous medium. In this context, a Neumann-series method for a heterogeneous medium has been developed, but assumes vacuum boundary conditions [34]. Using the approaches outlined in this manuscript, the Neumann-series method for heterogeneous medium could be extended to model reflection of photons at the boundary of the tissue. Further, in several bio-photonics imaging applications, reflection can occur between different tissues. For example, in transcranial imaging, refractive index changes occur between skull and cerebrospinal fluid in the brain. Improving the Neumann-series approach to model photon propagation in such media is another important area of future study. Solving a set of coupled equations using a treatment similar to Lehtikangas et al. [3] offers a mechanism for this purpose.
The reconstruction method proposed in this manuscript was used to estimate the absorption coefficient for a medium that had uniform optical coefficients. This was because our main purpose was to investigate the effect of modeling boundary conditions on the accuracy of estimating the optical coefficients. However, the proposed reconstruction method is general and an important direction of research is to refine the method for estimating the scattering and attenuation coefficients for a non-uniform medium. We are pursuing these efforts as part of an NIH BRAIN Initiative project on imaging in vivo neurotransmitter modulation of brain network activity in real time.
In this manuscript, we have focused on the deterministic approaches to solve the RTE. Another widely prevalent set of techniques for modeling photon transport are the stochastic MC-based techniques. The output of MC-based techniques is considered as the gold standard in several computational validation studies of the RTE, including in this manuscript. However, for image reconstruction, the use of stochastic MC-based techniques can lead to noisy estimates of the derivate or the need for simulating a large number of photons, which can be computationally expensive. In this context, several techniques have been proposed to accelerate the MC method [44]. However, using MC method to reconstruct the scattering coefficients is still challenging [45][46][47]. In contrast, deterministic approaches such as the Neumann-series RTE yield an analytical expression for the gradient that is not effected by noise. Given the trade-offs between the stochastic and deterministic RTE approaches, research on integrating these approaches for improved image reconstruction is an important frontier.

Conclusions
This work has proposed and investigated a Neumann-series-based radiative transport equation (RTE) for improved modeling of photon propagation in tissue. The method models photon reflection at the interface of tissue and external medium using Fresnel's equations. Further, the method was used to develop an algorithm to reconstruct the absorption and scattering coefficients of a scattering medium. Computational studies demonstrated that the method yielded more accurate modeling of photon transport for a 3D diffuse optical imaging (DOI) system in comparison to when the photon reflection was not modeled. In addition, the method yielded substantially more accurate estimates of the absorption coefficients for the 3D DOI system. The results demonstrate the importance of accounting for the reflection of photons at boundary when modeling photon propagation.
The above expression has a simple physical interpretation, namely that the photons are reflected back into the medium along the negativeẑ direction, as would be expected. Finally, applying the attenuation operator (Eq. (4)) yields (41)