We propose a formalism to incorporate boundary conditions in a Neumann-series-based radiative transport equation. The formalism accurately models the reflection 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 diffuse optical imaging. In computational studies, it was observed that the average root-mean-square error (RMSE) for the output images and the estimated absorption coefficients reduced by 38% and 84%, respectively, when the reflection boundary conditions were incorporated. These results demonstrate the importance of incorporating boundary conditions that model the reflection of photons at the tissue-external medium interface.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
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 , and 50% or more errors in estimating the optical coefficients of the tissue .
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–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 . 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–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–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 (PN) [19–22] require solving (N + 1)2 equations, where N is the number of Legendre polynomials. Similarly, the discrete ordinates methods (SN) [23, 24] requires solving N(N + 2) equations, where N is the number of angular directions. Limitations of the SN method also include false scattering and ray effects caused by spatial and angular discretization  and constraints on numerical parameters to ensure stability . Typically the vacuum boundary conditions were assumed in the PN methods [21, 22, 27] and early implementations of the SN methods [24, 28, 29]. More recently boundary conditions that model the reflection of photons were implemented for the SN methods [23, 30]. To overcome the computational complexity of the SN and PN methods, a simplified spherical harmonics (SPN) method has been proposed [10, 14, 31, 32]. The SPN methods implement partial-reflection boundary conditions. However, the SPN solutions also have several limitations .
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  and heterogeneous scattering media . 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 . 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].
2.1. 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 . In terms of photons, the photon distribution function 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 cm. Let denote a monochromatic mono-energetic emission source, and let denote the scattering phase function. For this source, the RTE can be written as 
One choice for scattering phase function in biological tissue is the Henyey-Greenstein function.
SimilarlyFig. 1. Using these notations, the integral form of the RTE is given by 
The equation has a very intuitive physical interpretation; the contribution of photons that have scattered n times is given by the term 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. , 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 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 is defined using Fresnel’s laws of reflection. Let , and denote the angle of incidence, reflection, and normal to the boundary surface. In accordance with Fresnel’s laws of reflection,
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:Fig. 3. For example, the term 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 . The term denotes photons that are scattered in the medium and subsequently reflected from a boundary surface.
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 . 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 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 
For this source term, the expression for the term is derived in Appendix A and is
The second term 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 . 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:
Denote the scattering operator in spherical harmonics by and the source term Ξ′ by χ′. The term is denoted in the spherical harmonic basis as . Using Eq. 16 yields
The third term 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 , where the index i denotes the different planes, denotes the unit normal vector to the ith plane and pi denotes the distance of the ith plane from the origin along the normal vector. With this notation, the expression for the term, as derived in Appendix B, is
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 mth pixel with center at rm by Φm. Assuming that the distribution function over a pixel is approximately constant, we get Appendix B.
2.3. 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:
Denote the mean noiseless image as a function of the scattering and absorption coefficient by . In our reconstruction approach, the objective is to estimate the value of µ that minimizes the norm of the error between the measured data g and the . Mathematically, the objective is to estimate that satisfies
The calculation of the above equation requires computing the derivative of with respect to µa and µs. Let denote the detector response function of the mth detector pixel. Then, the mth component of , denoted by , is
From Eq. 7,
As derived previously , terms containing partial derivatives of and are
Further, as derived in Appendix C
From Eq. (7), . Thus
Define . 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 .
From the gradient of the distribution function and using Eq. (27), the gradient of gm 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 remains the same as the above derivation.
2.4. 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 MCX Monte-Carlo (MC) technique . 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 gm,RTE and gm,MC denote the measurements acquired by the mth 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 . 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 Neumann-series-based reconstruction method that did not contain the reflection operator . 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 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  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%.
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 R2 << 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 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 , where i denotes the imaginary unit. Thus, the Neumann-series RTE can be solved for each frequency component by simply replacing µtot by . 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 . 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–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 µsL is around 5 , 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  and MC-DA method  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 . 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 . 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.  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 . However, using MC method to reconstruct the scattering coefficients is still challenging [45–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.
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.
Appendix A: derivation of the term
The photons described by the term are incident on the boundary defined by the plane z = L. We define this infinitesimally thin boundary by the expression δ(z − L). Applying the reflection operator, as defined by Eq. 8, on the distribution function within this boundary yields
Using the sifting property of the delta function and the fact that the normal to the boundary surface is the direction vector yields
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
Appendix B: derivation of the term
To derive the term, first note that, as derived previously 
Using Eq. 4, we obtain
Next, the reflection operator must be applied. This operator acts only on locations within a thin boundary region, denoted by . Applying the reflection operator, as defined by Eq. 8, on with this definition for the boundary yields
To simplify this integral, replace by r′, so that , λ |r − r′| and d3r′ = λ2dλdΩ. For notational simplicity, define
The effect of the attenuation operator (Eq. 4) on this term is
Using the scaling property of the delta function
For notational simplicity, defineEq. (48) in Eq. (47) and using the above definition for λ21i yields
Substituting the above expression in Eq. (22) yields the following expression for the transmitted flux detected by the m th pixel due to the XRXKXΞ term:
Appendix C: derivation of the term
Taking the derivative of the above expression with respect to µn yields
Separating the exponential integral into two parts, we obtain
NIH BRAIN Initiative Award R24 MH106083.
Dean Wong acknowledges contract work with Lilly, Lundbeck, Roche and Dart pharmaceuticals.
The authors thank Matthew Kupinski, PhD, Eric Clarkson, PhD and Harrison Barrett, PhD at the Center for Gamma Ray Imaging, University of Arizona, for helpful discussions.
References and links
1. R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]
2. G. W. Faris, “Diffusion equation boundary conditions for the interface between turbid media: a comment,” J. Opt. Soc. Amer. A 19, 519–520 (2002). [CrossRef]
3. O. Lehtikangas, T. Tarvainen, A. Kim, and S. Arridge, “Finite element approximation of the radiative transport equation in a medium with piece-wise constant refractive index,” J. Comput. Phys. 282, 345–359 (2015). [CrossRef]
4. M. Machida, G. Y. Panasyuk, J. C. Schotland, and V. A. Markel, “The Green’s function for the radiative transport equation in the slab geometry,” J. Phys. A: Math. and Theoretical 43, 065402 (2010). [CrossRef]
5. A. Liemert and A. Kienle, “Analytical approach for solving the radiative transfer equation in two-dimensional layered media,” J. Quant. Spectro. Radiative Transfer 113, 559–564 (2012). [CrossRef]
7. A. Pulkkinen and T. Tarvainen, “Truncated Fourier-series approximation of the time-domain radiative transfer equation using finite elements,” J. Opt. Soc. Am. A 30, 470–478 (2013). [CrossRef]
8. A. Gibson and H. Dehghani, “Diffuse optical imaging,” Phil. Tran. A. Math. Phys. Eng. Sci. 367, 3055–3072 (2009). [CrossRef]
9. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, 1–43 (2005). [CrossRef]
10. A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comp. Phys. 220, 441–470 (2006). [CrossRef]
11. A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. 43, 1285–1302 (1998). [CrossRef] [PubMed]
12. E. D. Aydin, C. R. de Oliveira, and A. J. Goddard, “A comparison between transport and diffusion calculations using a finite element-spherical harmonics radiation transport method,” Med. Phys. 29, 2013–2023 (2002). [CrossRef] [PubMed]
14. M. Chu, K. Vishwanath, A. D. Klose, and H. Dehghani, “Light transport in biological tissue using three-dimensional frequency-domain simplified spherical harmonics equations,” Phys. Med. Biol. 54, 2493–2509 (2009). [CrossRef] [PubMed]
15. R. Aronson, “Boundary conditions for diffusion of light,” J. Opt. Soc. Am. A 12, 2532–2539 (1995). [CrossRef]
16. A. Kienle and M. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A 14, 246–254 (1997). [CrossRef]
17. A. D. Kim, “Correcting the diffusion approximation at the boundary,” J. Opt. Soc. Am. A 28, 1007–1015 (2011). [CrossRef]
18. O. Lehtikangas, T. Tarvainen, and A. D. Kim, “Modeling boundary measurements of scattered light using the corrected diffusion approximation,” Biomed. Opt. Express 3, 552–571 (2012). [CrossRef] [PubMed]
19. J. K. Fletcher, “A solution of the neutron transport equation using spherical harmonics,” J. Phys. A: Math. Gen. 16, 2827–2835 (1983). [CrossRef]
20. K. Kobayashi, H. Oigawa, and H. Yamagata, “The spherical harmonics method for the multigroup transport equation in x-y geometry,” Ann. Nucl. Energy 13, 663–678 (1986). [CrossRef]
21. S. Wright, M. Schweiger, and S. Arridge, “Reconstruction in optical tomography using the PN approximations,” Meas. Sc. Tech. 18, 79 (2006). [CrossRef]
22. P. S. Mohan, T. Tarvainen, M. Schweiger, A. Pulkkinen, and S. R. Arridge, “Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements,” J. Comp. Phys. 230, 7364–7383 (2011). [CrossRef]
23. A. D. Klose, U. Netz, J. Beuthan, and A. H. Hielscher, “Optical tomography using the time-independent equation of radiative transfer - part 1: forward model,” J. Quant. Spectroscopy Rad. Transfer 72, 691–713 (2002). [CrossRef]
24. K. Ren, G. Bal, and A. H. Hielscher, “Frequency domain optical tomography based on the equation of radiative transfer,” SIAM J. Sc. Comp. 28, 1463–1489 (2006). [CrossRef]
25. P. Coelho, “The role of ray effects and false scattering on the accuracy of the standard and modified discrete ordinates methods,” J. Quant. Spect. Rad. Transfer 73, 231–238 (2002). [CrossRef]
26. M. Roger, C. Caliot, N. Crouseilles, and P. Coelho, “A hybrid transport-diffusion model for radiative transfer in absorbing and scattering media,” J. Comp. Phys. 275, 346–362 (2014). [CrossRef]
27. E. Aydin, C. De Oliveira, and A. Goddard, “A finite element-spherical harmonics radiation transport model for photon migration in turbid media,” J. Quant. Spectroscopy Rad. Transfer 84, 247–260 (2004). [CrossRef]
28. O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Problems 14, 1107 (1998). [CrossRef]
29. F. Asllanaj and S. Fumeron, “Applying a new computational method for biological tissue optics based on the time-dependent two-dimensional radiative transfer equation,” J. Biomed. Opt . 17, 075007 (2012). [PubMed]
30. F. Asllanaj, S. Contassot-Vivier, A. Liemert, and A. Kienle, “Radiative transfer equation for predicting light propagation in biological media: comparison of a modified finite volume method, the Monte Carlo technique, and an exact analytical solution,” J. Biomed. Opt . 19, 015002 (2014). [CrossRef]
31. Z. Yuan, X.-H. Hu, and H. Jiang, “A higher order diffusion model for three-dimensional photon migration and image reconstruction in optical tomography,” Phys. Med. Biol . 54, 65 (2009). [CrossRef]
32. Y. Lu, B. Zhu, H. Shen, J. C. Rasmussen, G. Wang, and E. M. Sevick-Muraca, “A parallel adaptive finite element simplified spherical harmonics approximation solver for frequency domain fluorescence molecular imaging,” Phys. Med. Biol . 55, 4625 (2010). [CrossRef] [PubMed]
33. A. K. Jha, M. A. Kupinski, T. Masumura, E. Clarkson, A. A. Maslov, and H. H. Barrett, “Simulating photon-transport in uniform media using the radiative transfer equation: A study using the neumann-series approach,” J. Opt. Soc. Am. A 29, 1741–1757 (2012). [CrossRef]
34. A. K. Jha, M. A. Kupinski, H. H. Barrett, E. Clarkson, and J. H. Hartman, “A three-dimensional Neumann-series approach to model light transport in nonuniform media,” J. Opt. Soc. Am. A 29, 1885–1899 (2012). [CrossRef]
35. A. K. Jha, Y. Zhu, J. Dreyer, J. Kang, A. Gjedde, D. Wong, and A. Rahmim, “Incorporating boundary conditions in the integral form of the radiative transfer equation for transcranial imaging,” in “Biomedical Optics 2016,” (2016), p. JW3A.47. [CrossRef]
36. A. K. Jha, Y. Zhu, D. F. Wong, and A. Rahmim, “A radiative transfer equation-based image-reconstruction method incorporating boundary conditions for diffuse optical imaging,” in SPIE Med. Imag. (2017).
37. A. K. Jha, “Retrieving Information from Scattered Photons in Medical Imaging,” Ph.D. thesis, College of Optical Sciences, University of Arizona, Tucson, AZ, USA (2013).
38. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22, 1779–1792 (1995). [CrossRef] [PubMed]
39. H. H. Barrett and K. J. Myers, Foundations of Image Science, 1st ed. (Wiley, 2004).
40. A. K. Jha, E. Clarkson, and M. A. Kupinski, “An ideal-observer framework to investigate signal detectability in diffuse optical imaging,” Biomed. Opt. Express 4, 2107–2123 (2013). [CrossRef] [PubMed]
42. T. Tarvainen, M. Vauhkonen, V. Kolehmainen, S. R. Arridge, and J. P. Kaipio, “Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and non-scattering regions,” Phys. Med. Biol. 50, 4913 (2005). [CrossRef] [PubMed]
43. L. Wang and S. L. Jacques, “Hybrid model of Monte Carlo simulation and diffusion theory for light reflectance by turbid media,” J. Opt. Soc. Am. A 10, 1746–1752 (1993). [CrossRef]
44. C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Optics 18, 050902 (2013). [CrossRef]
45. C. K. Hayakawa, J. Spanier, F. Bevilacqua, A. K. Dunn, J. S. You, B. J. Tromberg, and V. Venugopalan, “Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,” Opt. Lett. 26, 1335–1337 (2001). [CrossRef]
46. Y. P. Kumar and R. M. Vasu, “Reconstruction of optical properties of low-scattering tissue using derivative estimated through perturbation Monte-Carlo method,” J. Biomed. Optics 9, 1002–1012 (2004). [CrossRef]
47. R. Yao, X. Intes, and Q. Fang, “A rapid approach to build Jacobians for optical tomography via Monte Carlo method and photon “replay”,” in “Optics in the Life Sciences Congress,” (Optical Society of America, 2017), p. BoW3A.3. [CrossRef]