PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Biomed Opt. Author manuscript; available in PMC 2009 July 30.
Published in final edited form as:
J Biomed Opt. 2008; 13(2): 024016.
doi:  10.1117/1.2907168
PMCID: PMC2718535
NIHMSID: NIHMS122702

Integral equations of the photon fluence rate and flux based on a generalized Delta-Eddington phase function

Abstract

We present a generalized Delta-Eddington phase function to simplify the radiative transfer equation to integral equations with respect to both photon fluence rate and flux vector. The photon fluence rate and flux can be solved from the system of integral equations. By comparing to the Monte Carlo simulation results, the solutions of the system of integral equations accurately model the photon propagation in biological tissue over a wide range of optical parameters.

Keywords: photon propagation, radiative transfer equation, Delta-Eddington phase function, optical imaging

1 Introduction

The propagation of light through the biological tissue is a complicated process involving both absorption and scattering. The photon propagation model provides insight into the interaction between light and tissues and is essential for tomographic imaging with visible and near-IR light.1,2 The radiative transfer equation (RTE) is considered the golden standard for biomedical applications.1,2 Analytical solutions for the RTE are available for few simple geometries, and numerical solutions, such as the discrete ordinate method3,4 and the spherical harmonic method,5 often lead to enormous computational cost, especially to solve inverse problems such as optical tomography,6 bioluminescence tomography,710 and fluorescence tomography.11,12 Monte Carlo (MC) is a statistical simulation method in which the paths of photons are traced as they are scattered and absorbed within the medium.13 The MC method is well established to produce accurate estimates for light propagation in tissues. However, due to its statistical nature, the MC method has the disadvantage of requiring a long computation time; therefore, it is usually used as a reference method for other approaches. The popular diffusion approximation14,15 (DA) to RTE is widely used in the field of biophotonics because of its high computational efficiency. Nevertheless, DA assumes weakly anisotropic scattering and works well only in a highly scattering and weakly absorbing medium.16 It is also not suitable for modeling light propagation in the visible spectrum in which biological tissue present significant photon absorption.17

In the RTE, the phase functions describe the scattering characteristics of the medium and significantly influence the precision of the solution and the efficiency of the computation. Because the exact form of the phase function is currently unknown, the popular Henyey-Greenstein (HG) function18 and the Delta-Eddington function19 are usually used to approximate the true phase functions in biomedical applications. These two functions can be written in closed forms with a single free parameter g, called the anisotropic factor, which is often considered to be independent of the tissue scattering and absorption. Based on the Delta-Eddington phase function, a generalized diffusion model was presented to simulate photon propagation in highly absorbing medium and smaller source-detector separations.20 The inverse problem of optical parameters was presented based on the RTE with the Delta-Eddington function.21 The HG function was proven to be the most accurate in terms of the angular dependence of single scattering events in biological tissues.2,13 However, the RTE with the HG phase function is difficult to simplify further. In this paper, we present a generalized Delta-Eddington function to approximate the real phase function. Based on this new definition of Delta-Eddington phase function, the RTE can be reduced to a system of integral equations with respect to both the photon fluence rate and the flux vector. The solution of the system of integral equations enables a highly accurate prediction for the photon propagation in biological tissue over a wide range of optical parameters of biomedical interest.

2 Phase Approximation Model

The biological tissue scatters light strongly in the forward direction,22 and the scattering phase function can be modeled by a generalized Delta-Eddington function:19

p(υ·υ)=14π[(1f)(1+3hυ·υ)+2fδ(1υ·υ)],
(1)

where f [set membership] [−1, +1] is the weight factor measuring the anisotropy of the photon scattering, which we call the anisotropy weight; and h is a asymmetry factor of the phase function to modulate the weakly anisotropic scattering. The phase function is a linear combination of the weakly anisotropic scattering and the strongly peaked forward scattering. The original Delta-Eddington phase function rigidly defines the parameters f and h as g2 and g/(1+g), respectively.2,19 These parameters are related to a single free parameter g, where g is the anisotropic factor defined as the mean of the cosine of the scattering angles. Here, the proposed generalized Delta-Eddington function defines that the anisotropy weight and asymmetry factor are related to the photon absorbing and scattering coefficients in the medium. The optical properties of the medium are characterized by the anisotropy weight and asymmetry factor along with the absorption and scattering coefficients.

Based on the generalized Delta-Eddington phase function [Eq. (1)], the RTE

υ·L(r,υ)+(μa+μs)L(r,υ)=μsS2L(r,υ)p(υ,υ)dυ   +14πQ(r),   rΩ,
(2)

can be simplified as

υ·L(r,υ)+μ˜trL(r,υ)=μ˜s4π[Φ(r)+3hυ·J(r)]   +14πQ(r),   rΩ,
(3)

where Ω is the region of interest, L(r,υ) is the photon radiance (in watts mm−2 sr−1), Q(r) is the isotropic source (in watts mm−3), μs is the scattering coefficient (in mm−1), and μa is the absorption coefficient (in mm−1), [mu]s=(1−f)μs, and [mu]tra+[mu]s. The photon fluence rate Φ(r) and photon flux vector J(r) are, respectively, defined by

Φ(r)=S2L(r,υ)dυ   and   J(r)=S2υL(r,υ)dυ.
(4)

Equation (3) is a linear, first-order differential equation for the photon propagation in a heterogeneous medium, and the radiance L(r,υ) can be formulated as23

L(r,υ)=14π0R{μ˜s[Φ(rρυ)3hυ·J(rρv)]   +Q(rρυ)}exp[0ρμ˜tr(rtυ)dt]dρ   +L(rRυ,υ)exp[0Rμ˜tr(rtυ)dt],
(5)

where R is a scalar so that rRυ [set membership] [partial differential]Ω, representing the distance from point r to the boundary [partial differential]Ω along the direction υ. Integrating Eq. (5) over all the solid angles, we have

Φ(r)=14πS20R{μ˜s[Φ(rρυ)3hυ·J(rρυ)]   +Q(rρυ)}exp[0ρμ˜tr(rtυ)dt]dρdυ   +S2L(rRυ,υ)exp[0Rμ˜tr(rtυ)dt]dυ.
(6)

Multiplying both sides of Eq. (5) by the unit vector υ, and integrating Eq. (5) over all the solid angles, we obtain an equation for the flux vector:

J(r)=14πS20R{μ˜s[Φ(rρυ)3hυ·J(rρυ)]+Q(r   ρυ)}exp(0ρμ˜tr(rtυ)dt)υdρdυ+S2L(r   Rυ,υ)exp(0Rμ˜tr(rtυ)dt)υdυ.
(7)

In an optical imaging experiment, if no external photon enters the object Ω, L(rRυ,υ) in Eq. (6) and Eq (7) would vanish for the matched refractive indices on the boundary. However, the refractive index nissue in the biological tissues is higher than the refractive index nair of the surrounding air, and photons will be internally reflected at the boundary. In this case, L(rRυ,υ) describes the reflected radiance from a fraction of photons reflected on the boundary [partial differential]Ω, and can be expressed as

L(rRυ,υ)=rdL(rRυ,υ),   υ=υ2(υ·n)n,   υ·n<0,
(8)

where n is the outward unit normal at rRυ on [partial differential]Ω, and the internal reflection coefficient rd can be calculated14,16 by rd =−1.4399η−2+0.7099η−1+0.6681+0.0636η, η=ηissue/ ηair. From the first-order approximation of the radiance L(rRυ,υ′), the reflected radiance L(rRυ,υ) on the boundary can be approximated by

L(rRυ,υ)rd4π{Φ(rRυ)+3[υ2(υ·n)n]·J(r   Rυ)},  rRυΩ.
(9)

Substituting Eq. (9) into Eq. (6) and Eq (7) and performing a variable transformation from the polar coordinates r−ρυ to the Cartesian coordinates r′, we obtain a system of integral equations with respect to the photon fluence rate and flux vector:

{Φ(r)=14πΩ{μ˜s[Φ(r)3hυ·J(r)]+Q(r)}G(r,r)dr   +rd4πΩ[Φ(r)+3υ·J(r)]G(r,r)(υ·n)drJ(r)=14πΩ{μ˜s[Φ(r)3hυ·J(r)]+Q(r)}G(r,r)υdr   +rd4πΩ[Φ(r)+3υ·J(r)]G(r,r)(υ·n)υdr,
(10)

where

υ=rrrr,   υ=υ2(υ·n)n,   and   G(r,r)   =1|rr|2exp[0|rr|μ˜tr(rtυ)dt].

For simplicity, we call Eq. (10) the PA equation because it is derived from an approximate phase function. The PA equation is a well-posed system of integral equations of the second kind, and it enables an accurate prediction to the photon fluence rate and flux over a wide range of optical coefficients. The photon fluence rate presents the distribution of photon density in the spatial domain, while the photon flux describes the axial photon energy transfer and is directly related to experimental measurement data on the boundary. The exiting photon flux Γ(r) at r on the surface of object can be linked with the photon flux,2

Γ(r)=J(r)·n/(1rd),   for0<rd<1,   rΩ,
(11)

where [mu]s, [mu]t, and h (or, equivalently, μa, μs, f, and h) in the PA equation represent the optical parameters of the medium. When the absorption coefficient μa and scattering coefficient μs in the medium are unknown, we can directly determine the optical parameters [mu]s, [mu]t, and h using optical tomography techniques based on the PA equation.6 If the absorption coefficient μa, the scattering coefficient μs, and conventional anisotropic factor g in the tissue are known, based on a simple homogenous numerical phantom (such as a spherical or cylindrical object) with a known light source setting, the exiting photon flux on the boundary of the phantom can be generated using MC simulation with a appropriate phase function, such as the HG phase function. Then, an optimization procedure can be performed to determine the anisotropy weight f and asymmetry factor h by matching the exiting photon flux on the phantom boundary obtained from the PA equation to the MC-simulated counterpart. These exiting photon fluxes distribute on the boundary of the phantom and reflect the variation with different distances from the light source. Hence, the anisotropy weight f and asymmetry factor h are independent of the phantom geometry.

3 Numerical Experiments

To perform the numerical computation based on Eq. (10), the region of interest Ω can be discretized into finite elements with N vertex nodes, and the photon fluence rate Φ(r) and flux J(r) are approximated in terms of nodal-based basis functions24 ϕj(r)(j=1,2, … ,N),

{Φ(r)=j=1NΦ(rj)φj(r)J(r)=j=1NJ(rj)φj(r).
(12)

Substituting Eq. (12) into Eq. (10), we obtain a system of linear equations with respect to the photon fluence rate and flux:

{{Φ}=M1{Φ}+H1{J}+{Q1}{J}=M2{Φ}+H2{J}+{Q2},
(13)

where M1,M2,H1,H2, {Q1}, and {Q2} represent the corresponding discrete integral kernels and source vector in Eq. (10).

Extensive numerical experiments were conducted to validate the accuracy of the PA equation by MC simulation. The experiments are based on a cylindrical phantom of 20 mm diameter and 20 mm height. We set up a spherical source of 1.0 mm radius with a power of 10 nw at position of (−4.0,0.0,10.0) in the phantom. The phantom was then discretized into 25,335 tetrahedral elements and 4833 nodes, as shown in Fig. 1. A total of 1226 virtual detectors were allocated over the surface of the phantom to record the exiting photons. A wide range of optical parameters were respectively assigned to the phantom to mimic both high-absorption and low-absorption media. A reduced scattering albedo (RSA) defined by μs(1−g)/μa was used to measure the ratio of photon scatting to absorption. The relative refractive index η of the boundary was set at 1.37. An HG phase function with an anisotropic coefficient of 0.9 was assumed for the MC simulation. The MC simulation and the computational scheme of linear Eq. (13) were respectively performed to compute the photon fluence rate and exiting photon flux at detectors for three sets of optical parameters in Table 1. These results showed that the solutions of the PA equation were in excellent agreement with the results obtained from the MC simulation with a relative error below 3.8%. Here the relative error is defined as ΣdetectorsMC− ΓPA|/ΣdetectorsΓMC. Figure 2(a) to Figure 4(a) present a comparison between the PA and MC data for RSAs 3.57, 7.25, and 100, respectively. In contrast, we also performed the numerical experiments with the same phantom setting to compute the photon fluence rates based on the DA model. As a result, the relative errors between the photon fluence rate obtained from the DA model and from the MC simulation were as high as 31.2 and 14.6% for an RSAs of 3.57 and 7.25, respectively, as shown in Fig. 2(b) and Fig 3(b). As we had expected, the DA model had a satisfactory performance with relative errors of 3.8% for an RSA of 100, as shown in Fig. 4(b). The numerical experiments demonstrated that the solution obtained from the PA equation is more accurate than the results from the DA model as compared to the MC data over a broad range of optical parameters of biomedical interest. The computational time of the PA equation was about 45 min for both fluence rate and flux in our numerical experiments, while that of the MC simulation with the photon number of 107 was about 80 min, and that of the DA model in the same setting about 3 min. The running time is measured under a 2.8-GHz Intel Xeon CPU with 4 G bytes of memory.

Fig 1
Cylindrical phantom with a spherical light source discretized by tetrahedral elements.
Fig 2
Comparison between the MC simulation and the PA and DA models with the optical parameters μa=0.35 mm−1, μs=12.5 mm−1, h=0.9, f=0.938, and η =1.37. The detectors are sorted by the increasing order of the MC simulation ...
Fig 3
Comparison between the MC simulation and the PA and DA models with the optical parameters μa=0.20 mm−1, μs=14.5 mm−1, h=0.9, f=0.950, and η=1.37. The detectors are sorted by the increasing order of the MC simulation ...
Fig. 4
Comparison between the MC simulation and the PA and DA models with the optical parameters μa=0.008 mm−1, μs=8.0 mm−1, h=0.98, f=0.973, and η=1.37. The detectors are sorted by the increasing order of the MC simulation ...
Table 1
Optical parameters used in the numerical experiments.

4 Discussion and Conclusions

The phase function in the RTE describes the scattering characteristics of the medium and is strictly related to the scattering and absorption coefficients of the medium. The generalized Delta-Eddington phase function is characterized by anisotropy weight and asymmetry factor and shows that the anisotropy weight and asymmetry factor are related to the photon absorbing and scattering in the medium. The proposed concept improves the popular interpretations of the Delta-Eddington and HG phase functions, which are independent of absorbing and scattering in a medium. Based on the generalized Delta-Eddington phase function, the RTE can be significantly reduced to a system of integral equations, which enables the simultaneous computation of both the photon fluence rate and flux vector. MC simulation experiments demonstrated that the solution of the system of integral equations is highly accurate to model photon propagation over a wide range of optical parameters. The system of integral equations is suitable to model the photon propagation in heterogeneous media and can be applied to optical molecular imaging for a small animal with complex anatomical structures.

Acknowledgments

This work was partially supported by the National Institutes of Health under Grants EB001685, EB006036, EB008476, and CA127189.

References

1. Ishimaru A. Wave Propagation and Scattering in Random Media. Vol. 1. New York: Academic; 1978.
2. Welch AJ, van Gemert MJC. Optical and Thermal Response of Laser-Irradiated Tissue. New York: Plenum Press; 1995.
3. Chandrasekhar S. Radiative Transfer. New York: Dover Publications; 1960.
4. Abdoulaev GS, Hielscher AH. Three-dimensional optical tomography with the equation of radiative transfer. J. Electron. Imaging. 2003;12:594–601.
5. Klose AD, Larsen EW. Light transport in biological tissue based on the simplified spherical harmonics equations. J. Comput. Phys. 2006;220:441–470.
6. Arridge SR. Optical tomography in medical imaging. Inverse Probl. 1999;15:R41–R93.
7. Wang G, Hoffman EA, McLennan G, Wang LV, Suter M, Meinel J. Development of the first bioluminescent CT scanner. Radiology. 2003;229:566.
8. Wang G, Li Y, Jiang M. Uniqueness theorems in bioluminescence tomography. Med. Phys. 2004;31:2289–2299. [PubMed]
9. Cong W, Wang G, Kumar D, Liu Y, Jiang M, Wang LV, Hoffman EA, McLennan G, McCray PB, Zabner J, Cong A. Practical reconstruction method for bioluminescence tomography. Opt. Express. 2005;13:6756–6771. [PubMed]
10. Wang G, Cong W, Durairaj K, Qian X, Shen H, Sinn P, Hoffman E, McLennan G, Henry M. In vivo mouse studies with bioluminescence tomography. Opt. Express. 2006;14:7801–7809. [PubMed]
11. Zacharakis G, Ripoll J, Weissleder R, Ntziachristos V. Fluorescent protein tomography scanner for small animal imaging. IEEE Trans. Med. Imaging. 2005;24:878–885. [PubMed]
12. Joshi A, Bangerth W, Sevick-Muraca EM. Adaptive finite element based tomography for fluorescence optical imaging in tissue. Opt. Express. 2004;12:5402–5417. [PubMed]
13. Wang LH, Jacques SL, Zheng LQ. MCML—Monte Carlo modeling of photon transport in multi-layered tissues. Comput. Methods Programs Biomed. 1995;47:131–146. [PubMed]
14. Groenhuis RAJ, Ferwerda HA, Ten Bosch JJ. Scattering and absorption of turbid materials determined from reflection measurements. 1: Theory. Appl. Opt. 1983;22:2456–2467. [PubMed]
15. Keijzer M, Star WM, Storchi PRM. Optical diffusion in layered media. Appl. Opt. 1988;27:1820–1824. [PubMed]
16. Schweiger M, Arridge SR, Hiraoka M, Delpy DT. The finite element method for the propagation of light in scattering media: boundary and source conditions. Med. Phys. 1995;22:1779–1792. [PubMed]
17. Ripoll J, Yessayan D, Zacharakis G, Ntziachristos V. Experimental determination of photon propagation in highly absorbing and scattering media. J. Opt. Soc. Am. A. 2005;22:546–551. [PubMed]
18. Henyey L, Greenstein J. Diffuse radiation in the galaxy. Astrophys. J. 1941;93:70–83.
19. Joseph JH, Wiscombe WJ, Weinman JA. The Delta-Eddington approximation for radiative flux transfer. J. Atmos. Sci. 33;1976:2452–2459.
20. Venugopalan V, You JS, Tromberg BJ. Radiative transport in the diffusion approximation: an extension for highly absorbing media and small source-detector separations. Phys. Rev. E. 1998;58:2395–2407.
21. McCormick NJ. Inverse radiative transfer with a delta-Eddington phase function. Astrophys. Space Sci. 1987;129:331–334.
22. Cheong W-F, Prahl SA, Welch AJ. A review of the optical properties of biological tissues. IEEE J. Quantum Electron. 1990;26:2166–2185.
23. Natterer F, Wübbeling F. Mathematical Methods in Image Reconstruction. Philadelphia: SIAM; 2001.
24. Rao SS. The Finite Element Method in Engineering. Boston: Butterworth-Heinemann; 1999.