Translator Disclaimer
1 March 2008 Integral equations of the photon fluence rate and flux based on a generalized Delta-Eddington phase function
Author Affiliations +
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.

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, 7, 8, 9, 10 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

Eq. 1

p(vv)=14π[(1f)(1+3hvv)+2fδ(1vv)],
where f[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

Eq. 2

vL(r,v)+(μa+μs)L(r,v)=μsS2L(r,v)p(v,v)dv+14πQ(r),rΩ,
can be simplified as

Eq. 3

vL(r,v)+μ̃trL(r,v)=μ̃s4π[Φ(r)+3hvJ(r)]+14πQ(r),rΩ,
where Ω is the region of interest, L(r,v) is the photon radiance (in wattsmm2sr1 ), Q(r) is the isotropic source (in wattsmm3 ), μs is the scattering coefficient (in mm1 ), and μa is the absorption coefficient (in mm1 ), μ̃s=(1f)μs , and μ̃tr=μa+μ̃s . The photon fluence rate Φ(r) and photon flux vector J(r) are, respectively, defined by

Eq. 4

Φ(r)=S2L(r,v)dvandJ(r)=S2vL(r,v)dv.
Equation 3 is a linear, first-order differential equation for the photon propagation in a heterogeneous medium, and the radiance L(r,v) can be formulated as23

Eq. 5

L(r,v)=14π0R{μ̃s[Φ(rρv)3hvJ(rρv)]+Q(rρv)}exp[0ρμ̃tr(rtv)dt]dρ+L(rRv,v)exp[0Rμ̃tr(rtv)dt],
where R is a scalar so that rRvΩ , representing the distance from point r to the boundary Ω along the direction v . Integrating Eq. 5 over all the solid angles, we have

Eq. 6

Φ(r)=14πS20R{μ̃s[Φ(rρv)3hvJ(rρv)]+Q(rρv)}exp[0ρμ̃tr(rtv)dt]dρdv+S2L(rRv,v)exp[0Rμ̃tr(rtv)dt]dv.
Multiplying both sides of Eq. 5 by the unit vector v , and integrating Eq. 5 over all the solid angles, we obtain an equation for the flux vector:

Eq. 7

J(r)=14πS20R{μ̃s[Φ(rρv)3hvJ(rρv)]+Q(rρv)}exp(0ρμ̃tr(rtv)dt)vdρdv+S2L(rRv,v)exp(0Rμ̃tr(rtv)dt)vdv.
In an optical imaging experiment, if no external photon enters the object Ω , L(rRv,v) in Eqs. 6, 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(rRv,v) describes the reflected radiance from a fraction of photons reflected on the boundary Ω , and can be expressed as

Eq. 8

L(rRv,v)=rdL(rRv,v),v=v2(vn)n,vn<0,
where n is the outward unit normal at rRv on Ω , 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(rRv,v) , the reflected radiance L(rRv,v) on the boundary can be approximated by

Eq. 9

L(rRv,v)rd4π{Φ(rRv)+3[v2(vn)n]J(rRv)},rRvΩ.
Substituting Eq. 9 into Eqs. 6, 7 and performing a variable transformation from the polar coordinates rρv to the Cartesian coordinates r , we obtain a system of integral equations with respect to the photon fluence rate and flux vector:

Eq. 10

{Φ(r)=14πΩ{μ̃s[Φ(r)3hvJ(r)]+Q(r)}G(r,r)dr+rd4πΩ[Φ(r)+3vJ(r)]G(r,r)(vn)drJ(r)=14πΩ{μ̃s[Φ(r)3hvJ(r)]+Q(r)}G(r,r)vdr+rd4πΩ[Φ(r)+3vJ(r)]G(r,r)(vn)vdr,
where
v=rrrr,v=v2(vn)n,andG(r,r)=1rr2exp[0rrμ̃tr(rtv)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

Eq. 11

Γ(r)=J(r)n(1rd),for0<rd<1,rΩ,
where μ̃s , μ̃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 μ̃s , μ̃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) ,

Eq. 12

{Φ(r)=j=1NΦ(rj)φj(r)J(r)=j=1NJ(rj)φj(r).
Substituting Eq. 12 into Eq. 10, we obtain a system of linear equations with respect to the photon fluence rate and flux:

Eq. 13

{{Φ}=M1{Φ}+H1{J}+{Q1}{J}=M2{Φ}+H2{J}+{Q2},
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 20mm diameter and 20mm height. We set up a spherical source of 1.0mm radius with a power of 10nw 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(1g)μ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 detectorsΓMCΓPAdetectorsΓMC . Figures 2a to 4a 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 Figs. 2b and 3b . 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. 4b . 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 45min for both fluence rate and flux in our numerical experiments, while that of the MC simulation with the photon number of 107 was about 80min , and that of the DA model in the same setting about 3min . The running time is measured under a 2.8-GHz Intel Xeon CPU with 4G bytes of memory.

Fig. 1

Cylindrical phantom with a spherical light source discretized by tetrahedral elements.

024016_1_024802jbo1.jpg

Fig. 2

Comparison between the MC simulation and the PA and DA models with the optical parameters μa=0.35mm1 , μs=12.5mm1 , h=0.9 , f=0.938 , and η=1.37 . The detectors are sorted by the increasing order of the MC simulation photon fluence rates.

024016_1_024802jbo2.jpg

Fig. 3

Comparison between the MC simulation and the PA and DA models with the optical parameters μa=0.20mm1 , μs=14.5mm1 , h=0.9 , f=0.950 , and η=1.37 . The detectors are sorted by the increasing order of the MC simulation photon fluence rates.

024016_1_024802jbo3.jpg

Fig. 4

Comparison between the MC simulation and the PA and DA models with the optical parameters μa=0.008mm1 , μs=8.0mm1 , h=0.98 , f=0.973 , and η=1.37 . The detectors are sorted by the increasing order of the MC simulation photon fluence rates.

024016_1_024802jbo4.jpg

Table 1

Optical parameters used in the numerical experiments.

μa(mm−1) μs(mm−1) μs(1−g)∕μa g η f h
0.3512.503.570.91.370.9380.90
0.2014.507.250.91.370.9500.90
0.0088.01000.91.370.9730.98

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. 

A. Ishimaru, Wave Propagation and Scattering in Random Media, 1 Academic, New York (1978). Google Scholar

2. 

A. J. Welch and M. J. C. van Gemert, Optical and Thermal Response of Laser-Irradiated Tissue, (1995) Google Scholar

3. 

S. Chandrasekhar, “Radiative Transfer,” (1960) Google Scholar

4. 

G. S. Abdoulaev and A. H. Hielscher, “Three-dimensional optical tomography with the equation of radiative transfer,” J. Electron. Imaging, 12 594 –601 (2003). https://doi.org/10.1117/1.1587730 1017-9909 Google Scholar

5. 

A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys., 220 441 –470 (2006). https://doi.org/10.1016/j.jcp.2006.07.007 0021-9991 Google Scholar

6. 

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl., 15 R41 –R93 (1999). https://doi.org/10.1088/0266-5611/15/2/022 0266-5611 Google Scholar

7. 

G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. Meinel, “Development of the first bioluminescent CT scanner,” Radiology, 229 (P), 566 (2003). 0033-8419 Google Scholar

8. 

G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys., 31 2289 –2299 (2004). https://doi.org/10.1118/1.1766420 0094-2405 Google Scholar

9. 

W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. V. Wang, E. A. Hoffman, G. McLennan, P. B. McCray, J. Zabner, A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express, 13 6756 –6771 (2005). https://doi.org/10.1364/OPEX.13.006756 1094-4087 Google Scholar

10. 

G. Wang, W. Cong, K. Durairaj, X. Qian, H. Shen, P. Sinn, E. Hoffman, G. McLennan, and M. Henry, “In vivo mouse studies with bioluminescence tomography,” Opt. Express, 14 7801 –7809 (2006). https://doi.org/10.1364/OE.14.007801 1094-4087 Google Scholar

11. 

G. Zacharakis, J. Ripoll, R. Weissleder, and V. Ntziachristos, “Fluorescent protein tomography scanner for small animal imaging,” IEEE Trans. Med. Imaging, 24 878 –885 (2005). https://doi.org/10.1109/TMI.2004.843254 0278-0062 Google Scholar

12. 

A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express, 12 5402 –5417 (2004). https://doi.org/10.1364/OPEX.12.005402 1094-4087 Google Scholar

13. 

L. H. Wang, S. L. Jacques, and L. Q. Zheng, “MCML—Monte Carlo modeling of photon transport in multi-layered tissues,” Comput. Methods Programs Biomed., 47 131 –146 (1995). https://doi.org/10.1016/0169-2607(95)01640-F 0169-2607 Google Scholar

14. 

R. A. J. Groenhuis, H. A. Ferwerda, and J. J. Ten Bosch, “Scattering and absorption of turbid materials determined from reflection measurements. 1: Theory,” Appl. Opt., 22 2456 –2467 (1983). 0003-6935 Google Scholar

15. 

M. Keijzer, W. M. Star, and P. R. M. Storchi, “Optical diffusion in layered media,” Appl. Opt., 27 1820 –1824 (1988). 0003-6935 Google Scholar

16. 

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). https://doi.org/10.1118/1.597634 0094-2405 Google Scholar

17. 

J. Ripoll, D. Yessayan, G. Zacharakis, and V. Ntziachristos, “Experimental determination of photon propagation in highly absorbing and scattering media,” J. Opt. Soc. Am. A, 22 546 –551 (2005). https://doi.org/10.1364/JOSAA.22.000546 0740-3232 Google Scholar

18. 

L. Henyey and J. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J., 93 70 –83 (1941). https://doi.org/10.1086/144246 0004-637X Google Scholar

19. 

J. H. Joseph, W. J. Wiscombe, and J. A. Weinman, “The Delta-Eddington approximation for radiative flux transfer,” J. Atmos. Sci., 33 2452 –2459 (1976). https://doi.org/10.1175/1520-0469(1976)033<2452:TDEAFR>2.0.CO;2 0022-4928 Google Scholar

20. 

V. Venugopalan, J. S. You, and B. J. Tromberg, “Radiative transport in the diffusion approximation: an extension for highly absorbing media and small source-detector separations,” Phys. Rev. E, 58 2395 –2407 (1998). https://doi.org/10.1103/PhysRevE.58.2395 1063-651X Google Scholar

21. 

N. J. McCormick, “Inverse radiative transfer with a delta-Eddington phase function,” Astrophys. Space Sci., 129 331 –334 (1987). https://doi.org/10.1007/BF00642650 0004-640X Google Scholar

22. 

W.-F. Cheong, S. A. Prahl, and A. J. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron., 26 2166 –2185 (1990). https://doi.org/10.1109/3.64354 0018-9197 Google Scholar

23. 

F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction, (2001) Google Scholar

24. 

S. S. Rao, The Finite Element Method in Engineering, (1999) Google Scholar
©(2008) Society of Photo-Optical Instrumentation Engineers (SPIE)
Wenxiang Cong, Haiou Shen, Alex Cong, and Ge Wang "Integral equations of the photon fluence rate and flux based on a generalized Delta-Eddington phase function," Journal of Biomedical Optics 13(2), 024016 (1 March 2008). https://doi.org/10.1117/1.2907168
Published: 1 March 2008
JOURNAL ARTICLE
6 PAGES


SHARE
Advertisement
Advertisement
Back to Top