OSA's Digital Library

Virtual Journal for Biomedical Optics

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

  • Editor: Gregory W. Faris
  • Vol. 2, Iss. 6 — Jun. 13, 2007
« Show journal navigation

An improved Monte Carlo diffusion hybrid model for light reflectance by turbid media

Bin Luo and Sailing He  »View Author Affiliations


Optics Express, Vol. 15, Issue 10, pp. 5905-5918 (2007)
http://dx.doi.org/10.1364/OE.15.005905


View Full Text Article

Acrobat PDF (243 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

This paper introduces an improved diffusion model which is accurate and fast in computation for the cases of μa/μ’s < 007 as good as the conventional diffusion model for the cases of μa/μ’s < 0.007 for surface measurement, hence more suitable than the conventional model to be the forward model used in the image reconstruction in the diffuse optical tomography. Deviation of the diffusion approximation (DA) on the medium surface is first studied in the Monte Carlo (MC) diffusion hybrid model for reflectance setup. A modification of DA and an improved MC diffusion hybrid model based on this modified DA are introduced. Numerical tests show that for media with relatively strong absorption the present modified diffusion approach can reduce the surface deviation significantly in both the hybrid and pure diffusion model, and consumes nearly no more computation time than the conventional diffusion approach.

© 2007 Optical Society of America

1. Introduction

Diffuse optical tomography (DOT) [1–2

1. A. Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today 48, 34–40 (1995). [CrossRef]

] using near-infrared light can image the optical properties of biological tissues non-invasively, therefore is used and promising in many biomedical applications such as breast cancer detection [3

3. K. T. Moesta, S. Fantini, H. Jess, S. Totkas, M. A. Franceschini, M. Kaschke, and P. M. Schlag, “Constrast features of breast cancer in frequency-domain laser scanning mammography,” J. Biomed. Opt. 3, 129–136 (1998) [CrossRef]

] and functional brain imaging [4–5

4. A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. 20, 435–442 (1997). [CrossRef] [PubMed]

]. The image reconstruction in DOT based on the theory of photon transport or photon density wave [6

6. A. Ishimaru, Wave propagation and scattering in random media, (Academic Press, New York1978).

] is a non-linear ill-posed inverse problem which requires an accurate and fast (in computation) forward model [2

2. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. Dimarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process Mag. 57–75 (2001). [CrossRef]

]. The radiance transport equation (RTE) [6

6. A. Ishimaru, Wave propagation and scattering in random media, (Academic Press, New York1978).

] and Monte Carlo (MC) simulation [7–9

7. B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys. 10, 824–830 (1983). [CrossRef] [PubMed]

] are two accurate models for light propagation in turbid medium. However, both the two models are time-consuming hence not suitable to be the forward model for the DOT image reconstruction. Thus, for highly scattering media with relatively weak absorption, people use the diffusion approximation (DA) to simplify the RTE to a more efficient (on computation) model, i.e. the diffusion model represented by the diffusion equation (DE) [6

6. A. Ishimaru, Wave propagation and scattering in random media, (Academic Press, New York1978).

]. Due to the DA, the DE suffers from some limitations. The first limitation is that the medium must be scattering dominated (i.e., the ratio of the absorption coefficient over the reduced scattering coefficient must be much less than 1, e.g., when, μ a/μs ≤ 0.005, the deviation of the DE solution from the RTE solution is usually smaller than 1%). Another limitation of DE is that it can not correctly model light propagation in the “near-source region” (the region near the incident point of the collimated light source) [10–13

10. 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–2462 (1983). [CrossRef] [PubMed]

]. This deviation of DE in the near-source region is referred as “near-field deviation” in the present paper. To overcome this difficulty, hybrid models that use the diffusion model in the far-from-source region and an accurate model in the near-source region was proposed, e.g., Monte Carlo diffusion hybrid model (MC simulation is used in the near-source region) [14–16

14. 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]

] and transport diffusion hybrid model (the RTE is used in the near-source region) [17–18

17. G. Bal and Y. Maday, “Coupling of transport and diffusion models in linear transport theory,” Math. Model Numer. Anal. 36, 69–86 (2002). [CrossRef]

].

However, in biomedical applications, it is often met that μ a/μs < 0.007 is not satisfied. For example, in human brain near-infrared imaging, μ a/μs of the scalp is about 0.0095, and μ a/μs of the gray matter is about 0.0164 [5

5. E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. II,” Appl. Opt. 42, 2915–2922 (2003). [CrossRef] [PubMed]

]. For those even less scattering dominated cases, the conventional diffusion model (no matter pure or hybrid) is not accurate enough to be the forward model used in the DOT image reconstruction. As an example, consider a setup of light reflectance as shown in Fig. 1, when μ a/μs ≥ 0.05, the deviation of the diffusion modeled light intensity at the detector from the RTE modeled value is usually greater than 7%. Figure 2 shows a numerical example of the comparison between the conventional MC diffusion hybrid model and the pure MC model. From Fig. 2 one sees that the deviation of the surface diffuse intensity is greater than 7% when the distance from the incident point is larger than 1.2 cm.

Fig. 1. A schematic diagram of the photon transportation in a turbid medium and the detection of the diffuse light on the surface.
Fig. 2. The pure MC model and the conventional MC diffusion hybrid model are used to evaluate the diffuse intensity on the surface of a semi-infinite homogeneous medium, and the results by the two models are compared. This figure shows the relative deviation of the surface diffuse intensity by the conventional MC diffusion hybrid model (from the pure MC model). The geometry of the setup is shown in Fig. 1. The optical parameters of the medium are: n = 1.35, g = 0.8, μs = 49.4cm-1, μa = 0.6cm-1. In the hybrid model in this example, the surface diffuse intensity in the region ρ ≥ 0.45cm is evaluated with the finite element method based on the conventional diffusion equation.

In this paper, we introduce an improved diffusion model which is accurate and fast in computation for the cases of μ a/μs < 0.07 as good as the conventional diffusion model for the cases of μ a/μs < 0.007 for surface measurement, hence more suitable than the conventional model to be the forward model used in the DOT image reconstruction. The idea is presented in the following two paragraphs.

The validity of the DA relies on that the diffuse radiance in the medium be nearly isotropic. However, due to the strong anisotropy of the physical material at the interface between the medium and the air (e.g., the mismatch of the refractive indices), it is reasonable that the diffuse light near the surface of the medium is more anisotropic than the light deep inside the medium. Therefore, the deviation of the diffuse modeled light intensity from the accurate value is larger on the surface than in the deep-inside region [19–20

19. 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]

]. This relatively larger deviation of the diffusion model on the surface is called the “surface deviation” in the present paper. Note that the measurements in DOT are usually taken on the surface of the biological medium, the surface deviation is a disadvantage for the conventional diffusion model to be the forward model for DOT.

By modifying some basic equations in the diffusing theory, we obtain an improved diffusion model. This model reduces the deviation on the surface to the level of the deviation deep inside the medium, therefore, is more accurate than the conventional diffusion model.

The accuracy of the improved model for the cases of μ a/μs < 0.07 is as good as the conventional diffusion model for the cases of μ a/μs < 0.007. Meanwhile, this improved model (as a forward model) consumes nearly no more computation time than the conventional diffusion model. Therefore, this model is more suitable than the conventional model to be the forward model for the image reconstruction in DOT.

To modify the conventional diffusion model to a model that is suitable for the cases of μ a/μs < 0.07, we first study in Section 2 the accuracy behavior of the basic equations (listed in Section 2.1) of the conventional DA when the conventional model is applied to the cases of 0.05 ≤ μ a/μs < 0.07. This study shows which equations listed in Section 2.1 are more inaccurate and should be modified to correctly model the light propagation in the near-surface region. In Section 3.1 several empirical formulae are proposed to modify the basic equations of the conventional DA. These modified basic equations represent the “modified DA”. A modified MC diffusion hybrid model (for light reflectance by turbid media) based on the modified DA is introduced in Section 3.2 and numerical tests for the modified model are provided in Section 3.3. Supported by the numerical experiments including that shown in Section 3.3, we give our conclusion in Section 4.

2. Surface deviation of the conventional diffusion approximation

2.1. The conventional diffusion approximation

In the present paper, we assume the incident light is from a continuous-wave (CW) source and thus the RTE is time-independent. In the far-from-source region of a scattering dominated medium, the light propagation is completely diffusive and governed by the following diffusion equation (far-from-source DE) [6

6. A. Ishimaru, Wave propagation and scattering in random media, (Academic Press, New York1978).

]

[(13μt(r))+μa(r)]Ud(r)=0,rΩ
(1)

where μ’t = μ’s + μa is the reduced transport coefficient, Ud(r)14π4πJd(r,s)dω is the average light intensity, Jd(r,s) is the diffuse radiance (Wsr-1cm-2) at point r in the direction of unit vector s, and Ω represents the far-from-source region of the turbid medium. Note that Eq. (1) (called the “conventional DE” in the present paper) is derived from the following two equations [6

6. A. Ishimaru, Wave propagation and scattering in random media, (Academic Press, New York1978).

]:

μa(r)Ud(r)+14πFd(r)=0
(2)
14μμt(r)Fd(r)+13Ud(r)=0
(3)

where F d (r) ≜ ∫4π Jd(r, s)sdω is the diffuse flux.

On the surface of the turbid medium, the partial current boundary condition can be expressed as [21–22

21. M. Keijzer, W. M. Star, and P. R. M. Storchi, “Optical diffusion in layered media,” Appl. Opt. 27, 1820–1824 (1988). [CrossRef] [PubMed]

]

Ud(r)=A2πnFd(r),rS
(4)

where S represents the surface of the medium, n is the outward normal to S, and A is calculated with the following expression [22

22. D. Contini, F. Martelli, and G. Zaccanti, ”Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

]

A=1+3π0π2R(θ)cos2(θ)sin(θ)dθ120π2R(θ)cos(θ)sin(θ)dθ
(5)

where R(θ) is the Fresnel reflection coefficient for the diffuse light hitting the medium surface (from inside) at an angle of θ (with respect to n). Coefficient A is a function of the relative refractive index n of the medium and can be approximated with a polynomial fit of n in practical calculation [22

22. D. Contini, F. Martelli, and G. Zaccanti, ”Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

].

Equations (2)-(5) are all based on the following equation for an isotropic medium [6

6. A. Ishimaru, Wave propagation and scattering in random media, (Academic Press, New York1978).

]:

JdrsUd(r)+34πFd(r)s
(6)

However, the diffuse light is more anisotropic in the region near the surface (as compared with other regions of the medium), Therefore, Eq. (6) is relatively less accurate in the near-surface region, leading to the “surface deviation” of DA.

In the following Section 2.2, we introduce briefly three implementations of the conventional MC diffusion hybrid model as the tools to study the surface deviation.

2.2. Three implementations of the conventional MC diffusion hybrid model

Figure 3 shows a simple example of a Gaussian beam of light impinging normally on the surface of a semi-infinite turbid medium. The denotations V 1, V 2, ABCD, EFGH and ÉF́ǴH́ are specified in the caption of Fig. 3. Base on these five denotations, we define the following denotations:

{Ω=V1V2,volumeofV1withV2excavatedΓ0=V1ABCD,boundaryareaofV1exceptABCDΓ1=ABCDEFGH,areaofABCDwithEFGHexcavatedΓ2=EFGHEFGH,areaofEFGHwithEFGHexcavatedΓ3=ΩˉV2=V2EFGH,boundaryareaofV2exceptEFGH
Fig. 3. A thin beam (e.g., a Gaussian beam) of light incidents orthogonally onto the surface (the plane marked with “ABCD”) of a semi-infinite turbid medium. V 1 is a cubic volume (whose surface is square ABCD) of the medium. V 2 is a small cube (whose surface is square EFGH) included in V 1. Square E’F’G’H’ which includes EFGH is slightly larger than EFGH. Square ABCD, EFGH and E’F’G’H’ are all centered at the incident point O. MC simulation is used to evaluate the diffuse intensity in V 2 and on square E’F’G’H’. Finite element method is used to solve the diffusion equation in the volume of V 1 with the small cube V 2 excavated.

The denotations Γi (i = 0,1,2,3) are defined for the convenience of describing the various boundary conditions to be used in the present paper.

In the conventional MC diffusion hybrid model for the example shown in Fig. 3, we first use MC simulation to obtain the diffuse light distribution in the near-source region (including V 2 and Γ2). Then, the finite element method (FEM) [23–25

23. S. C. Brenner and L. R. Scott, Mathematical Theory of Finite Element Methods, Texts in Appl. Math. 15, (New York, Springer-Verlag1994).

] based on the conventional DE [Eq. (1)] is applied only in the far-from-source region Ω (The size of V 2 is chosen large enough so that the light in Ω is completely diffusive.) The diffuse intensity at point r in Ω decreases rapidly when the distance between r and O (the incident point) increases. Volume V 1 is chosen large enough so that the diffuse intensity on Γ0 is small enough and zero boundary condition can be applied at Γ0 for the DE solution.

Below we describe three implementations of the conventional MC diffusion hybrid model (called “Conventional Hybrid I”, “Conventional Hybrid II” and “Conventional Hybrid III”) using three different boundary conditions on Γ2. In all these three implementations, the FEM is used to solve Eq. (1) (combined with different boundary conditions).

The boundary conditions for Eq. (1) used in Conventional Hybrid I are

{Ud(r)=0,rΓ013μt(r)nUd(r)+12AUd(r)=0,rΓ1Γ213μt(r)nUd(r)=14πnFd(r),rΓ3
(7)

where the second equation in the above system is the so-called “partial current boundary condition” [21–22

21. M. Keijzer, W. M. Star, and P. R. M. Storchi, “Optical diffusion in layered media,” Appl. Opt. 27, 1820–1824 (1988). [CrossRef] [PubMed]

], which is derived from Eqs. (3) and (4). The third equation in system (7) is derived from Eq. (3), and term 14πnFd(rΓ3) can be obtained from the MC simulation result in the near-source region.

In Conventional Hybrids II and III, we use the same boundary conditions for Eq. (1) on Γ0, Γ1 and Γ3 as in Conventional Hybrid I. For the boundary condition on Γ2, we use in Conventional Hybrid II

13μt(r)nUd(r)=14πnFd(r),rΓ2
(8)

and use in Conventional Hybrid III

Ud(r)=α(r),rΓ2
(9)

where α(r ∈ Γ2) can be obtained from the MC simulation. One way to evaluate α(r ∈ Γ2) with MC simulation is to evaluate Jd(r ∈ Γ2,s) at first, then calculate α(r ∈ Γ2) with the definition of Ud:α(rΓ2)=14π4πJd(rΓ2,s)dω. Actually, this integration can be evaluated by the MC simulation itself.

The advantage of the Dirichlet constraint (9) is that it guarantees accurate value α(r ∈ Γ2) for Ud(r ∈ Γ2), and hence improves the accuracy of Ud(r) for r near Γ2. Note that the square ÉF́ǴH́ does not need to be much bigger than square EFGH, and square EFGH should be set a proper size. In the MC diffusion hybrid models, Ud on Γ2 is evaluated with MC simulation. The computation time for the simulation to evaluate Ud on Γ2 under a required precision increases exponentially with the distance between Γ2 and the source. Empirically, this distance is set about 5/μt́, i.e., the size of the “near-source region” is about 10/μt́. Outside this near-source region, the near-field deviation is quite small and it is suitable to use the diffusion model instead of the MC model to save the computation time. Setting a larger size (more than 10μt́) of the near-source region will consume extra time for the computation. For Conventional Hybrids I and III, one can set ÉF́ǴH́ almost identical to EFGH (i.e., Γ2 reduces to the four edge lines of square EFGH).

For ordinary usage, Conventional Hybrids I-III are applied to the cases of μa/μś < 0.007 . But in the following subsection, we apply Conventional Hybrids I-III to a case of μa/μś = 0.06 and compare their solutions to study the behavior of the conventional DA in the cases of μa/μś ≥ 0.05 .

2.3. Numerical study for the surface deviation

Here we study the surface deviation with a numerical example using the setup shown in Fig. 3. The medium is assumed to be homogeneous semi-infinite with the same optical parameters as listed in the caption of Fig. 2. Square EFGH (upper surface of V 2, Fig. 3) and ÉF́ǴH́ are set by

{squareEFGH={(x,y,z)x0.4380cm,y0.4380cm,z=0},squareEFGH={(x,y,z)x0.4517cm,y0.4517cm,z=0}

In this example, the diffuse intensity on the medium surface is evaluated by four different approaches: the pure MC simulation, Conventional Hybrids I, II and III. Figure 4 shows the numerical results for the deviation (from the result of the pure MC simulation) of the diffuse intensity on the positive x-axis obtained by Conventional Hybrids I-III. The relative deviations are defined by

DEVHybridI,II,III(x)=sfUdHybridI,II,III(x)sfUdMC(x)1
(10)

where sfUd MC(x), sfUd Hybrid I, II, III(x) are the surface diffuse intensities at position (x,0,0) obtained by the pure MC simulation and Conventional Hybrid (I, II or III), respectively.

Fig. 4. The relative deviations (from the result of the pure MC simulation) of the diffuse intensity on the positive x-axis obtained by Conventional Hybrids I-III

Note that the only difference among Conventional Hybrids I-III is the boundary condition used on Γ2. Therefore, as shown in Fig. 4, the significant difference between DEV hybrid II(x) and DEV hybrid III(x) in region 0.4380cm ≤ x ≤ 0.4517cm (i.e., sfUd hybrid II(r ∈ Γ2) is much greater than sfUd MC(r ∈ Γ2) while sfUd hybrid III(r ∈ Γ2) is close to sfUd MC(r ∈ Γ2)) implies that Eq. (8) is not quite accurate, i.e., Eq. (3) for r near the surface is not quite accurate in this example. Also note that DEV hybrid II(x) is a bit worse than DEV hybrid I(x). This indicates that Eq. (8) is less accurate than the second equation in system (7) [derived from Eqs. (3) and (4)]. This difference between DEV hybrid II(x) and DEV hybrid I(x) implies that Eq. (4) is not accurate. The underlying reason for the inaccuracy of Eqs. (3) and (4) is the relatively strong anisotropy of the diffuse light near the surface of the medium.

From the above analysis one sees that one way to reduce the surface deviation is to modify Eqs. (3) and (4).

3. The modified DA and the improved diffusion (hybrid or pure) models

3.1. The modified diffusion approximation

In the present paper we propose to modify the conventional DA [described by Eq. (2)-Eq. (4)] by rewriting Eq. (3) and Eq. (4) into the following two equations:

14πμt(r)Fd(r)+1+η1(r)3Ud(r)=0
(3’)
Ud(r)=[1+η2(r)]A2πnFd(r),rS
(4’)

where η 1 (r) and η 2 (r) are modification parameter functions introduced to reduce the surface deviation. They are usually selected empirically. For the setup and the coordinate system shown in Fig. 3 (the source beam incidents along z axis, and x-y plane is defined on the surface of the medium), we propose that η 1(r) and η 2(r) be approximated as:

{η1(r)c1(n)μa(r)μt(r)exp[0zμt(x,y,z)dz],η2(r)c2(n)μa(r)μt(r),rS
(11)

where c 1 (n) and c 2 (n) are constants determined only by the refractive index n of the medium.

Here we give an explanation for the use of the η 1(r) and η 2(r) in Eq. (11). 1. For light reflectance setup shown in Fig. 3, the deviation of the surface measurement evaluated by the conventional diffusion model from that by the MC model increases almost proportionally with μa/μś [22

22. D. Contini, F. Martelli, and G. Zaccanti, ”Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

]. This correlativity is also verified by the numerical experiments done by us. Therefore, the modification should be direct proportion to μa/μś. 2. The anisotropic light becomes nearly isotropic when it travels a path length in the scattering dominated medium. Therefore, the modification should decrease rapidly with the depth from the surface. It’s natural to choose η 1(r) to decrease exponentially with the dimensionless path length. 3. The deviation of the DE on the surface increases with the mismatch of the refractive index at the interface [19-20

19. 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]

]. Thus c 1 (n) and c 2 (n) should be determined by the refractive index n. For a fixed n, c 1 and c 2 are two constants. Usually their values are unknown, however they are not arbitrary. The values of c 1 and c 2 for a fixed n can be estimated from large numbers of numerical examples. In this paper, we use n = 1.35. We have done hundreds of numerical comparisons between the modified DA and the pure MC simulation for various optical properties of turbid media. In each numerical experiment, we choose various pairs of values of c 1 and c 2 to find a pair of values which makes the modified DA consistent with the pure MC simulation best. We found that for the case of n = 1.35, when we use the values c 1 = 3.14 and c 2 = 2.2, the modified DA is consistent with the pure MC simulation very good in most cases. Though the best pair values of c 1 and c 2 are unknown, they should be near the pair values of c 1 = 3.14 and c 2 = 2.2 which are used for the present modified DA in this paper.

3.2. The improved MC diffusion hybrid model (the present hybrid model

The diffusion equation and the boundary condition used in the proposed MC diffusion hybrid model (called “Present Hybrid model” hereafter) are modified as

[(D(r))+μa(r)]Ud(r)=0,rΩ
(12)

and

{Ud(r)=0,rΓ0D(r)nUd(r)+12[1+η2(r)]AUd(r)=0,rΓ1Ud(r)=α(r),rΓ2D(r)nUd(r)=14πnFd(r),rΓ3
(13)

where α(r ∈ Γ2) and 14πnFd(rΓ3) are obtained with the MC simulation, and D(r) is the diffusion coefficient that is defined in the present paper by

D(r)=1+η1(r)3μt(r)
(14)

Equation (12) is derived from Eqs. (2) and (3’). The second equation in system (13) (called the “modified partial current boundary condition”) is derived from Eqs. (3’) and (4’). Note that when η 1 (r) = 0, The definition of D(r) by Eq. (14) reduces to the conventional definition of the diffusion coefficient D(r) = [3μt́(r)]-1 (Ref. [6

6. A. Ishimaru, Wave propagation and scattering in random media, (Academic Press, New York1978).

]) and at the same time the modified DE [Eq. (12)] reduces to the conventional DE [Eq. (1)]. Absorption-independent diffusion coefficient D(r) = [3μt́(r)]-1 (corresponding to η 1(r) = μa(r)/μś(r)) was suggested in Refs. [22

22. D. Contini, F. Martelli, and G. Zaccanti, ”Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

] and [26

26. T. Durduran, A. G. Yodh, B. Chance, and D. A. Boas, “Does the photon-diffusion coefficient depend on absorption,” J. Opt. Soc. Am. A 14, 3358–3365 (1997). [CrossRef]

]. We modify D(r) by using the η 1(r) function described in Eq. (11). The partial current boundary condition in the second equation of system (13) is modified by the η 2(r) in Eq. (11). And in the present paper, we use the values of c 1 = 3.14 and c 2 = 2.2 , i.e., for light reflectance setups shown in Figs. 1 and 3, we propose the following diffusion coefficient and partial current boundary condition: (The refractive index of the medium n = 1.35)

{D(r)={1+3.14μa(r)μt(r)exp[0z(x,y,z)dz]}[3μt(r)]Ud(r)=12π[1+2.2μa(r)μt(r)]AnFd(r),rS
(15)

The diffusion approximation modified by Eq. (15) is referred as the “present modified DA” hereafter in the present paper. Our numerical experiments have verified that the present modified DA is more accurate than the DA modified with D(r) = [3μś(r)]-1 in Refs. [22

22. D. Contini, F. Martelli, and G. Zaccanti, ”Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

] and [26

26. T. Durduran, A. G. Yodh, B. Chance, and D. A. Boas, “Does the photon-diffusion coefficient depend on absorption,” J. Opt. Soc. Am. A 14, 3358–3365 (1997). [CrossRef]

]. We use an example in Section 3.3 to illustrate this (Fig. 7).

In the FEM, Eqs. (12) and (13

13. T. Spott and L. O. Svaasand, “Collimated light sources in the diffusion approximation,” Appl. Opt. 39, 6453–6465 (2000). [CrossRef]

) are usually transformed into the following Galerkin variation equation for Ud(r):

Ω[D(r)ϕ(r)Ud(r)+μa(r)ϕ(r)Ud(r)]dV+Γ112[1+η2(r)]Aϕ(r)Ud(r)dσ=Γ3ϕ(r)14π(n).Fd(r)dσ
(16)

where ϕ(r) is the test function satisfying ϕ(r ∈ Γ0 ∪ Γ2) = 0 , and Ud (r) is imposed with Ud(r ∈ Γ0) = 0 and Ud(r ∈ Γ2) = α(r ∈ Γ2).

A standard procedure of the FEM [23-25

23. S. C. Brenner and L. R. Scott, Mathematical Theory of Finite Element Methods, Texts in Appl. Math. 15, (New York, Springer-Verlag1994).

] can then be implemented to discretize Eq. (16) into a linear system of equations, which can be solved by an iterative algorithm such as the conjugated gradient iteration.

Note that when η 1 (r) ≡ 0 and η 2 (r) = 0 , Eq. (16) reduces to the Galerkin variation equation in the conventional model. This implies that the computation complexity of the forward problem in the modified model is the same as in the conventional model. Therefore, the modified DA approach consumes nearly no more computation time than the conventional DA approach.

3.3. Numerical tests for the modified DAs

{D(r)=[3μs(r)]1Ud(r)=12πAnFd(r),rS
(17)

This modification (Eq. (17), called “Modified DA II” in the present paper) corresponds to η 1 (r) = μa (r)/μś (r) and η 2 (r) = 0. A modification between Modified DA I and II is

{D(r)=[3μs(r)]1Ud(r)=12π[1+2.2μa(r)μt(r)]AnFd(r),rS
(18)

This modification (Eq. (18), called “Modified DA III” in the present paper) corresponds to η 1 (r)= μa(r)/μś(r) and η 2(r) = 2.2μa(r)/μt́(r).

Obviously, the modified DA can be applied in not only the hybrid model but also the pure diffusion model. In the pure diffusion model, Eq. (12) (the far-from-source DE) is replaced by the whole domain DE:
[(D(r))+μa(r)]Ud(r)=S(r),rV1
(19)
where V1 = V2 ∪ Ω is the whole domain of the medium, S(r) is the diffuse source which is derived from the reduced incident radiance [6

6. A. Ishimaru, Wave propagation and scattering in random media, (Academic Press, New York1978).

, 10

10. 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–2462 (1983). [CrossRef] [PubMed]

, 21

21. M. Keijzer, W. M. Star, and P. R. M. Storchi, “Optical diffusion in layered media,” Appl. Opt. 27, 1820–1824 (1988). [CrossRef] [PubMed]

]. In the far-from-source region Ω, S(r) is very small and can be considered as zero, i.e. S(r ∈ Ω) ≈ 0. And the boundary conditions system (13) is replaced by
{Ud(r)=0,rΓ0D(r)nUd(r)+12[1+η2(r)]AUd(r)=0,rsquareABCD
(20)

To show the accuracy of the present modified DA (Modified DA I), we give several numerical examples of the comparisons among various diffusion models named in Table 1.

Table 1. Nomenclature (used in Figs. 4–7) for the diffusion models

table-icon
View This Table

In all the examples the same geometry of the setup and the Cartesian coordinate system (shown in Figs. 1 and 3) are used. The source beam incidents along z axis, and x-y plane is defined on the surface of the medium. The incident point O is at (0, 0, 0).

The first four examples (Fig. 5) are used to show the advantage of the Present Hybrid model over the Conventional Hybrid I-III models. In each of the four examples, the medium is homogeneous with various optical parameters, and five different numerical approaches (the pure MC simulation, Conventional Hybrids I, II, III, and Present Hybrid) are used to evaluate the diffuse intensity distribution on the surface of the medium. The outputs of Conventional Hybrids I-III and Present Hybrid are compared with the output of the MC simulation. As shown in Fig. 5, the proposed hybrid model (Present Hybrid) is much more accurate than Conventional Hybrid I, II and III. In the first three examples [Figs. 5(a), 5(b) and 5(c)], Present Hybrid improves the accuracy significantly. In the fourth example [Fig. 5 (d)], Conventional Hybrid III is almost as accurate as Present Hybrid because of the small value of μa/μt́ (0.005) used in this example.

Fig. 5. Four MC diffusion hybrid models (Conventional Hybrids I-III and Present Hybrid) are used to evaluate the diffuse intensity on the surface of semi-infinite homogeneous media, and their results are compared with the pure MC simulation. The deviations are shown on the x axis in the region {(x, 0, 0) ∈ Γ1 ∪ Γ2}. Comparisons are done for four sets of medium parameters: (a) n = 1.35, g = 0.8, μs = 49.4cm-1, μa = 0.6cm-1; (b) n = 1.35, g = 0.9, μs = 99.4cm-1, μa = 0.6cm-1; (c) n = 1.35, g = 0.7, μs = 49.1cm-1, μa = 0.9cm-1; (d) n = 1.35, g = 0.8, μt = 99.9cm-1, μa = 0.1cm-1.

We use the fifth numerical example to show the accuracy improvement by the present modified DA in both the pure diffusion and hybrid models. In this example, five models are used to evaluate the diffuse intensity distribution on the surface of a semi-infinite inhomogeneous medium: the pure MC model, Conventional Pure diffusion model, Present Pure diffusion model (based on the present modified DA), Conventional Hybrid III, and Present Hybrid. In this example, the FEM is applied in the two pure diffusion models and in the far-from-source region for the two hybrid models. The medium is inhomogeneous and cylindrically symmetric about the z axis (the axis of the incident beam) and described by

{n1.35,g0.8,μs49.4cm1,μa(x,y,z)={[0.6+0.2cos((r0.4)π)]cm1,r<2.0cm0.4cm1,r2.0cm

where r = √x 2 + y 2 + z 2 ,(z ≥ 0). As shown in Fig. 6, Present Hybrid model is the most accurate (the deviation of Present Hybrid from the pure MC model is the smallest). The present pure diffusion model (denoted by “Present Pure” in Fig. 6) is much more accurate than the conventional pure diffusion model (denoted by “Conventional Pure” in Fig. 6), and in the far-from-source region, Present Pure is also more accurate than Conventional Hybrid III. The large values of DEV Conventional, Present Pure(x ≥ 0.1cm) are due to the near-field deviation of the diffusion model. Figure 6 shows that the present modified DA improves the accuracy significantly not only for the hybrid model but also for the pure diffusion model.

Fig. 6. Five models (Conventional Pure, Present Pure, Conventional Hybrid III, Present Hybrid, and the pure MC) are used to calculate the diffuse intensity on the surface of a semi-infinite inhomogeneous medium with cylindrical symmetry, and the results of the first four models are compared with the fifth model (the pure MC simulation). In the two hybrid models, the surface diffuse distribution in the region x ≤ 0.4cm is evaluated with MC simulation.
Fig. 7. (a). The deviations of the modified DAs in the pure diffusion model; (b). The deviations of the modified DAs in the MC diffusion hybrid model. The medium parameters are the same as used for Fig. 2, and the waist radius of the incident Gaussian beam is 0.03cm.

The last two examples (Fig. 7) are used to show the advantage of the present modified DA (Modified DA I) over Modified DA II (suggested in Refs. [22

22. D. Contini, F. Martelli, and G. Zaccanti, ”Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

] and [26

26. T. Durduran, A. G. Yodh, B. Chance, and D. A. Boas, “Does the photon-diffusion coefficient depend on absorption,” J. Opt. Soc. Am. A 14, 3358–3365 (1997). [CrossRef]

]) and Modified DA III. For the same medium as used in Figs. 2, 4 and 5(a), we apply Modified DA I, II and III in the pure diffusion model and the MC diffusion hybrid model described by Eqs. (12)-(14

14. 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]

) to evaluate the diffuse intensity on the surface of the medium. The numerical results are compared with the pure MC simulation. As shown in Fig. 7, the present modified DA (Modified DA I) is more accurate than Modified DA II and III in both the pure diffusion model (Fig. 7(a), the pure diffusion model applied with the present modified DA is marked with “Present Pure”) and the MC diffusion hybrid model (Fig. 7(b), the hybrid model applied with the present modified DA is marked with “Present Hybrid”).

As a summary of this subsection, Figures 5 and 6 demonstrate the significant improvement of the present modified DA [described by Eq. (15)] compared with the conventional DA; Fig. 5(d) is used to show that the modified DA reduces to the conventional DA when μa/μt́ → 0; and Fig. 7 shows that the present modified DA is also better than other modified DAs reported before (e.g., Modified DA II which was suggested in Refs. [22

22. D. Contini, F. Martelli, and G. Zaccanti, ”Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

] and [26

26. T. Durduran, A. G. Yodh, B. Chance, and D. A. Boas, “Does the photon-diffusion coefficient depend on absorption,” J. Opt. Soc. Am. A 14, 3358–3365 (1997). [CrossRef]

]).

At last, it should be pointed out that, the computation time for the FEM to solve the modified DA problems is almost the same long as for the FEM to solve the conventional DA problems.

4. Conclusion

For media with relatively strong absorption, the traditional diffusion approximation (DA) is not so accurate near the medium surface. Neither a pure nor hybrid conventional DA approach can reduce this deviation on far-from-source surface. We have introduced a method with some empirical formulae to modify the conventional DA (including the partial current boundary condition) in order to reduce the surface deviation. Numerical tests have shown that the improved models (pure diffusion or MC diffusion hybrid) which employ the present modified DA are much more accurate than the conventional models. With the combined consideration that the modified DA approaches consume nearly no more computation time than the conventional DA approaches, we conclude that the present improved models are more efficient hence more suitable than the conventional models to be forward models used in the image reconstruction in DOT.

Acknowledgments

This work was partially supported by the Science Foundation of Zhejiang Province, China (grant No. Y104255) and the National Science Foundation of China (grant No. 60372032). One of the authors (Bin Luo) would like to thank Dr. Jun Li for his valuable discussions.

References and links

1.

A. Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today 48, 34–40 (1995). [CrossRef]

2.

D. A. Boas, D. H. Brooks, E. L. Miller, C. A. Dimarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process Mag. 57–75 (2001). [CrossRef]

3.

K. T. Moesta, S. Fantini, H. Jess, S. Totkas, M. A. Franceschini, M. Kaschke, and P. M. Schlag, “Constrast features of breast cancer in frequency-domain laser scanning mammography,” J. Biomed. Opt. 3, 129–136 (1998) [CrossRef]

4.

A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. 20, 435–442 (1997). [CrossRef] [PubMed]

5.

E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. II,” Appl. Opt. 42, 2915–2922 (2003). [CrossRef] [PubMed]

6.

A. Ishimaru, Wave propagation and scattering in random media, (Academic Press, New York1978).

7.

B. C. Wilson and G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys. 10, 824–830 (1983). [CrossRef] [PubMed]

8.

L. Wang, S. L. Jacques, and L. Zhen, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47,131–146 (1995). [CrossRef] [PubMed]

9.

D. A. Boas, J. P. Culver, J. J. Stott, and A. K. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10,159–170 (2002). [PubMed]

10.

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–2462 (1983). [CrossRef] [PubMed]

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.

S. R. Arridge, “Optical tomography in medical imaging,” Inv. Probl. 15, R41–R93 (1999). [CrossRef]

13.

T. Spott and L. O. Svaasand, “Collimated light sources in the diffusion approximation,” Appl. Opt. 39, 6453–6465 (2000). [CrossRef]

14.

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]

15.

G. Alexandrakis, T. J. Farrrell, and M. S. Patterson, “Monte Carlo diffusion hybrid model for photon migration in a two-layer turbid medium in the frequency domain,” Appl. Opt. 39, 2235–2244 (2000). [CrossRef]

16.

T. Hayashi, Y. Kashio, and E. Okada, “Hybrid Monte Carlo-diffusion method for light propagation in tissue with a low-scattering region,” Appl. Opt. 42, 2888–2896 (2003). [CrossRef] [PubMed]

17.

G. Bal and Y. Maday, “Coupling of transport and diffusion models in linear transport theory,” Math. Model Numer. Anal. 36, 69–86 (2002). [CrossRef]

18.

T. Tarvainen, M. Vauhkonen, V. Kolehmainen, and J. Kaipio, “Finite element model for the coupled radiative transfer equation and diffusion approximation,” Int. J. Numer. Methods Eng. , 65, 383–405 (2006). [CrossRef]

19.

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]

20.

B. Q. Chen, K. Stamnes, and J. J. Stamnes, “Validity of the diffusion approximation in bio-optical imaging”, Appl. Opt. 40, 6356–6366 (2001). [CrossRef]

21.

M. Keijzer, W. M. Star, and P. R. M. Storchi, “Optical diffusion in layered media,” Appl. Opt. 27, 1820–1824 (1988). [CrossRef] [PubMed]

22.

D. Contini, F. Martelli, and G. Zaccanti, ”Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. 36, 4587–4599 (1997). [CrossRef] [PubMed]

23.

S. C. Brenner and L. R. Scott, Mathematical Theory of Finite Element Methods, Texts in Appl. Math. 15, (New York, Springer-Verlag1994).

24.

P. G. Cialet, Finite Element Method for Elliptic Problems, (North-Holland1978).

25.

J. M. Jin, Finite Element Method in Electromagnetics, 2nd ed. (New York, John Wiley & Sons2002).

26.

T. Durduran, A. G. Yodh, B. Chance, and D. A. Boas, “Does the photon-diffusion coefficient depend on absorption,” J. Opt. Soc. Am. A 14, 3358–3365 (1997). [CrossRef]

OCIS Codes
(170.5270) Medical optics and biotechnology : Photon density waves
(170.6960) Medical optics and biotechnology : Tomography
(170.7050) Medical optics and biotechnology : Turbid media
(290.1990) Scattering : Diffusion

ToC Category:
Scattering

History
Original Manuscript: January 8, 2007
Revised Manuscript: March 14, 2007
Manuscript Accepted: April 25, 2007
Published: April 30, 2007

Virtual Issues
Vol. 2, Iss. 6 Virtual Journal for Biomedical Optics

Citation
Bin Luo and Sailing He, "An improved Monte Carlo diffusion hybrid model for light reflectance by turbid media," Opt. Express 15, 5905-5918 (2007)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-15-10-5905


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. A. Yodh and B. Chance, "Spectroscopy and imaging with diffusing light," Phys. Today 48,34-40 (1995). [CrossRef]
  2. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. Dimarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, "Imaging the body with diffuse optical tomography," IEEE Signal Process Mag.57-75 (2001). [CrossRef]
  3. K. T. Moesta, S. Fantini, H. Jess, S. Totkas, M. A. Franceschini, M. Kaschke, and P. M. Schlag, "Constrast features of breast cancer in frequency-domain laser scanning mammography," J. Biomed. Opt. 3,129-136 (1998) [CrossRef]
  4. A. Villringer and B. Chance, "Non-invasive optical spectroscopy and imaging of human brain function," Trends Neurosci. 20,435-442 (1997). [CrossRef] [PubMed]
  5. E. Okada and D. T. Delpy, "Near-infrared light propagation in an adult head model. II," Appl. Opt. 42, 2915-2922 (2003). [CrossRef] [PubMed]
  6. A. Ishimaru, Wave propagation and scattering in random media, (Academic Press, New York 1978).
  7. B. C. Wilson and G. Adam, "A Monte Carlo model for the absorption and flux distributions of light in tissue," Med. Phys. 10, 824-830 (1983). [CrossRef] [PubMed]
  8. L. Wang, S. L. Jacques, and L. Zhen, "MCML - Monte Carlo modeling of light transport in multi-layered tissues," Comput. Methods Programs Biomed. 47,131-146 (1995). [CrossRef] [PubMed]
  9. D. A. Boas, J. P. Culver, J. J. Stott, and A. K. Dunn, "Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head," Opt. Express 10,159-170 (2002). [PubMed]
  10. 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-2462 (1983). [CrossRef] [PubMed]
  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. S. R. Arridge, "Optical tomography in medical imaging," Inv. Probl. 15,R41-R93 (1999). [CrossRef]
  13. T. Spott and L. O. Svaasand, "Collimated light sources in the diffusion approximation," Appl. Opt. 39,6453-6465 (2000). [CrossRef]
  14. 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]
  15. G. Alexandrakis, T. J. Farrrell, and M. S. Patterson, "Monte Carlo diffusion hybrid model for photon migration in a two-layer turbid medium in the frequency domain," Appl. Opt. 39,2235-2244 (2000). [CrossRef]
  16. T. Hayashi, Y. Kashio, and E. Okada, "Hybrid Monte Carlo-diffusion method for light propagation in tissue with a low-scattering region," Appl. Opt. 42,2888-2896 (2003). [CrossRef] [PubMed]
  17. G. Bal and Y. Maday, "Coupling of transport and diffusion models in linear transport theory," Math. Model Numer. Anal. 36,69-86 (2002). [CrossRef]
  18. T. Tarvainen, M. Vauhkonen, V. Kolehmainen, and J. Kaipio, "Finite element model for the coupled radiative transfer equation and diffusion approximation," Int. J. Numer. Methods Eng., 65,383-405 (2006). [CrossRef]
  19. 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]
  20. B. Q. Chen, K. Stamnes, and J. J. Stamnes, "Validity of the diffusion approximation in bio-optical imaging", Appl. Opt. 40,6356-6366 (2001). [CrossRef]
  21. M. Keijzer, W. M. Star, and P. R. M. Storchi, "Optical diffusion in layered media," Appl. Opt. 27,1820-1824 (1988). [CrossRef] [PubMed]
  22. D. Contini, F. Martelli, and G. Zaccanti, "Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory," Appl. Opt. 36,4587-4599 (1997). [CrossRef] [PubMed]
  23. S. C. Brenner and L. R. Scott, Mathematical Theory of Finite Element Methods, Texts in Appl. Math. 15, (New York, Springer-Verlag 1994).
  24. P. G. Cialet, Finite Element Method for Elliptic Problems, (North-Holland 1978).
  25. J. M. Jin, Finite Element Method in Electromagnetics, 2nd ed. (New York, John Wiley & Sons 2002).
  26. T. Durduran, A. G. Yodh, B. Chance, and D. A. Boas, "Does the photon-diffusion coefficient depend on absorption," J. Opt. Soc. Am. A 14,3358-3365 (1997). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited