1. 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–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–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 (P_N) [19–22] require solving (N + 1)² equations, where N is the number of Legendre polynomials. Similarly, the discrete ordinates methods (S_N) [23, 24] requires solving N(N + 2) equations, where N is the number of angular directions. Limitations of the S_N method also include false scattering and ray effects caused by spatial and angular discretization [25] and constraints on numerical parameters to ensure stability [26]. Typically the vacuum boundary conditions were assumed in the P_N methods [21, 22, 27] and early implementations of the S_N methods [24, 28, 29]. 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].

2. Methods

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 $\widehat{\mathbf{s}}$ by $w\left(\mathbf{r},\widehat{\mathbf{s}}\right)$. In terms of photons, the photon distribution function $w\left(\mathbf{r},\widehat{\mathbf{s}}\right)\mathrm{\Delta }V\mathrm{\Delta }\mathrm{\Omega }$ 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 $\widehat{\mathbf{s}}$ 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 $\mathrm{\Xi }\left(\mathbf{r},\widehat{\mathbf{s}}\right)$ denote a monochromatic mono-energetic emission source, and let $f\left(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime };\mathbf{r}\right)$ denote the scattering phase function. For this source, the RTE can be written as [33]

One choice for scattering phase function in biological tissue is the Henyey-Greenstein function.

Similarly

 

Fig. 1 Physical interpretation of the various operators in the RTE in considered DOI setup

Download Full Size | PDF

The equation has a very intuitive physical interpretation; the contribution of photons that have scattered n times is given by the term $\mathcal{X}{\left(\mathcal{K}\mathcal{X}\right)}^n\mathrm{\Xi }$ 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 $ℛ$ 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:

 

Fig. 2 A schematic describing the thin boundary layer in the considered DOI setup

Download Full Size | PDF

The operator $ℛ$ is defined using Fresnel’s laws of reflection. Let ${\widehat{\mathbf{s}}}^{\prime },\widehat{\mathbf{s}}$, and $\widehat{\mathbf{n}}$ 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 Physical interpretation of the three reflection terms in the RTE

Download Full Size | PDF

2.2. 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 $\mathcal{X}ℛ\mathcal{X}\mathrm{\Xi }$ 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 $\mathcal{X}ℛ\mathcal{X}\mathrm{\Xi }$ term is derived in Appendix A and is

The second term $\mathcal{X}\mathcal{K}\mathcal{X}ℛ\mathcal{X}\mathrm{\Xi }$ 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:

Denote the term $\mathcal{X}ℛ\mathcal{X}\mathrm{\Xi }$, 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 $\mathcal{D}$ and the source term Ξ′ by χ′. The term $\left[\mathcal{K}{\mathrm{\Xi }}^{\prime }\right]\left(\mathbf{r},\widehat{\mathbf{s}}\right)$ is denoted in the spherical harmonic basis as ${\left[\mathcal{D}{\chi }^{\prime }\right]}_{lm}\left(\mathbf{r}\right)$. 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 $\mathcal{X}ℛ\mathcal{X}\mathcal{K}\mathcal{X}\mathrm{\Xi }$ 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 ${\displaystyle {\sum }_i\delta \left(\mathbf{r}\cdot {\widehat{\mathbf{n}}}_i-p_i\right)}$, where the index i denotes the different planes, ${\widehat{\mathbf{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 $\mathcal{X}ℛ\mathcal{X}\mathcal{K}\mathcal{X}\mathrm{\Xi }$ 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 m^th pixel with center at r_m by Φ_m. Assuming that the distribution function over a pixel is approximately constant, we get [33]

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 $\overline{\mathbf{g}}\left(\mathit{\mu }\right)$. In our reconstruction approach, the objective is to estimate the value of µ that minimizes the ${\mathbb{L}}_2$ norm of the error between the measured data g and the $\overline{\mathbf{g}}\left(\mathit{\mu }\right)$. Mathematically, the objective is to estimate $\widehat{\mathit{\mu }}$ that satisfies

The calculation of the above equation requires computing the derivative of $\overline{\mathbf{g}}\left(\mathit{\mu }\right)$ with respect to µ_a and µ_s. Let $h_m\left(\mathbf{r},\widehat{\mathbf{s}}\right)$ denote the detector response function of the m^th detector pixel. Then, the m^th component of $\overline{\mathbf{g}}\left(\mathit{\mu }\right)$, denoted by ${\overline{g}}_m$, is

From Eq. 7,

As derived previously [40], terms containing partial derivatives of $\mathcal{X}$ and $\mathcal{K}$ are

Further, as derived in Appendix C

Substituting the expressions from Eqs. (29)–(33) in Eq. (28) and rearranging the terms yields

From Eq. (7), $\mathcal{X}\mathrm{\Xi }+\mathcal{X}ℛw+\mathcal{X}\mathcal{K}w=w$. Thus

Define $\mathcal{S}=-c_m{\varphi }_{n}w+\mathit{ϵ}{\varphi }_n{\mathcal{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 $\mathcal{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 $\frac{\partial {\overline{g}}_m}{\partial \mu }$ 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⁻¹ and 0.01 cm⁻¹. 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.

 

Fig. 4 The experimental setup simulated to evaluate the proposed framework

Download Full Size | PDF

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 [41]. 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,RTE 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 Neumann-series-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.

3. 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%.

 

Fig. 5 Comparison of the linear profiles of the images obtained with the MC and proposed Neumann-series RTE formalism, when the refractive index of the scattering medium is (a) 1.2 (b) 1.3 (c) 1.4 and (d) 1.5.

Download Full Size | PDF

 

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.

Download Full Size | PDF

Table 1. RMSE in the output image obtained using the Neumann-series RTE method with and without the reflection boundary conditions

View Table

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⁻¹ 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. 7 (a) Plot of estimated absorption coefficient as a function of the refractive index using the proposed reconstruction method and another method that does not model photon reflection. Error bars are also shown. (b) Plot of the RMSE between the true and estimated absorption coefficient as a function of the refractive index using the two reconstruction methods

Download Full Size | PDF

4. 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² << 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.

 

Fig. 8 The Fresnel reflection coefficient values as a function of the incident angle for various values of refractive index n, when the surrounding media has a refractive index of 1

Download Full Size | PDF

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 $\mathrm{\Xi }\left(\mathbf{r},\widehat{\mathbf{s}}\right)$ 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 ${\mu }_{\text{tot}}+\frac{i\nu }{\text{cm}}$, where i denotes the imaginary unit. Thus, the Neumann-series RTE can be solved for each frequency component by simply replacing µ_tot by ${\mu }_{\text{tot}}+\frac{i\nu }{\text{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–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 [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–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.

5. 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.