Photon diffusion near the point-of-entry in anisotropically scattering turbid media

doi:10.1038/ncomms1599

Nat Commun. Author manuscript; available in PMC 2012 December 13.

Published in final edited form as:

Nat Commun. 2011 December 13; 2: 587.

doi: 10.1038/ncomms1599

PMCID: PMC3370954

NIHMSID: NIHMS352343

Photon diffusion near the point-of-entry in anisotropically scattering turbid media

Edward Vitkin,^‡ Vladimir Turzhitsky,^‡ Le Qiu,^‡ Lianyu Guo, Irving Itzkan, Eugene B. Hanlon, and Lev T. Perelman^*

Center for Advanced Biomedical Imaging and Photonics Department of Obstetrics, Gynecology and Reproductive Biology and Department of Medicine Beth Israel Deaconess Medical Center Biological and Biomedical Sciences Program Harvard University Boston, Massachusetts 02215 USA

^‡These authors contributed equally to this work.

^*Correspondence and requests for materials should be addressed to L.T.P. (Email: ltperel/at/caregroup.harvard.edu).

The publisher's final edited version of this article is available at Nat Commun

See other articles in PMC that cite the published article.

Abstract

From astronomy to cell biology, the manner in which light propagates in turbid media has been of central importance for many decades. However, light propagation near the point-of-entry (POE) in turbid media has never been analytically described, until now. Here we report a straightforward and accurate method that overcomes this longstanding, unsolved problem in radiative transport. Our theory properly treats anisotropic photon scattering events and takes the specific form of the phase function into account. As a result, our method correctly predicts the spatially dependent diffuse reflectance of light near the POE for any arbitrary phase function. We demonstrate that the theory is in excellent agreement with both experimental results and Monte Carlo simulations for several commonly used phase functions.

Recent advances in optical spectroscopy of human epithelium, high resolution imaging and early cancer detection demand an accurate description of light propagation near the point where light enters the tissue^1–8. However, an analytical method capable of describing the propagation of photons near their point of entry (POE) in turbid media has eluded researchers in diverse fields for generations^9–12. Propagation of electrons in semiconductors¹³, anisotropic light diffusion in liquid crystals¹⁴, and coherent backscattering¹⁵ are only a few examples that can benefit from such a method. The most accurate treatment of such problems, when interference can be neglected, uses the radiative transport equation,

(1)

where L(r, s) is the radiance of light at position r traveling in the direction of unit vector s, S(r, s) describes the distribution of internal sources, and μ_a and μ_s are the absorption and scattering coefficients of the medium, respectively. The unity normalized function p(s, s′) is the phase function of the medium, which describes the angular distribution for a single scattering event (here, s′ is the incoming and s the outgoing direction of scatter)⁹.

Unfortunately, no analytical solution to this equation has been discovered for practical scenarios. An often used numerical method is Monte Carlo simulation. Monte Carlo simulation gives accurate results but does not provide any physical insight and often requires excessive computational time. Recent GPU based Monte Carlo methods significantly reduce computational time^16–19 thus extending the utility of Monte Carlo simulations to practical applications, but inverse problems with a large number of parameters will greatly benefit from the analytical approach provided here. In the absence of suitable numerical methods, the most widely used approach is to obtain an approximate solution using the diffusion approximation (DA) to the radiative transport equation. The DA works best when scattering is isotropic or near isotropic²⁰. It is not sensitive to the particular form of the phase function, reducing phase function information to a single parameter, equation M2

, called the average cosine of the phase function or the anisotropy factor, which in turn is incorporated into the reduced scattering coefficient, equation M3

. Hence,

becomes the only parameter in the DA to describe scattering. While the diffusion approximation can provide analytical solutions in some cases, it does not provide accurate results for the critically important region near the POE.

This limitation of the DA has been a major impediment, particularly for light scattering in biological media. Many investigators have attempted to amend the DA by solving the radiative transport equation using more complicated phase functions^21–23. One such approach, called the δ-P₁ approximation, is to solve the equation for the delta-Eddington phase function, which results in some improvements for characterizing forward scattering media²⁴. Another approach is to decompose the phase function into a series of Legendre polynomials and to construct a system of radiative transport equations corresponding to each Legendre term. Because it becomes very cumbersome to solve the system of partial differential equations as the number of Legendre terms increases, the P₃ approximation uses only the 0^th through 3^rd terms in the Legendre series²⁵. Both approximations are substantial improvements over the DA because they take the phase function into account and broaden the parameter space with which the approximation can be used^26,27. However, until now, no proposed analytical method has addressed the important regime of light scattering near the POE (e.g., at source-detector separations less than equation M5

) well enough to be widely used as common practice.

Here, we introduce a novel approach that results in the first successful attempt to analytically describe spatially dependent diffuse reflectance near the POE in turbid media with arbitrary phase function and absorption. The phase function corrected diffusion approximation (PFC) takes advantage of the delta-isotropic phase function, described below, to derive a two-part expression for the diffuse reflectance. The first part is identical to the DA for a pencil beam illumination, while the second part approximates the changes in reflectance due to the specific form of the phase function. The resulting predictions are shown to be in excellent agreement with experiment and also with Monte Carlo simulations at all distances, including the critical area near POE.

Results

Phase function corrected photon diffusion

We consider the problem of calculating the reflectance of light from a semi-infinite turbid medium with no internal sources (S(r, s = 0) illuminated with a narrow collimated beam (see Fig. 1), an important, general problem often treated using the DA.

Figure 1

Semi-infinite turbid medium geometry

The radiance in this problem can be divided into three main parts: an on-axis part due to the narrow collimated beam, a predominately isotropic, phase function independent part and an anisotropic, phase function dependent part. Similarly, any arbitrary phase function can be expressed as a sum of a delta function, an isotropic part and an anisotropic part. We achieve this by expressing an arbitrary phase function, p(s· s′), as the sum of the delta-isotropic phase function with the same anisotropy factor and a deviation term

(2)

where p_DI(s · s′)=[1− g +2g ·†(1−s·s′)]/4π is the delta-isotropic phase function²⁸. This decomposition of the phase function enables us to separate the radiance into the three analogous parts described above.

By substituting equation (2) into equation (1) and using the reduced scattering coefficient, we can express equation (1) as

(3)

Without the second integral on the right-hand side, equation (3) would be the equation of radiative transport for the case of an isotropically scattering medium, which, as mentioned above, can be described using the DA. Thus, indeed, it is logical to present the radiance L(r, s) as the sum of a predominantly isotropic diffuse part, L_d(r, s) a phase function dependent part, (L_p(r, s) and an on-axis part due to the narrow collimated beam, L_c(r, s).

After substitution for L(r, s) in equation (3), we seek an expression for L_c(r, s) which satisfies the boundary condition for a narrow collimated beam and a portion of equation (3), equation M8

. Note that this portion is described by equation M9

rather than μ_s and specifically accounts for low-angle scattered photons²⁹ (see Fig. 1 for illustration). The boundary condition corresponding to a narrow collimated beam entering a semi-infinite medium is equation M10

for inward directions ( equation M11

), where p₀ is the power of the incoming beam, ρ and z are the cylindrical coordinates and equation M12

is a unit vector along the z axis. This yields³⁰ equation M13

and the boundary conditions equation M14

and

for

The transport equation (3) can then be written as a system of two equations for the remaining parts of the radiance L_d(r, s) and L_p(r, s). It is natural then to distribute the two remaining L_c(r, s) dependent terms so that the isotropic, phase function independent term equation M17

is combined is combined with the equation for the predominately isotropic L_d(r, s) and the anisotropic, phase function dependent term equation M18

is combined with the equation for the anisotropic L_p(r, s).

The first equation,

(4)

leads to the standard diffusion approximation and the second equation,

(5)

retains all of the phase function dependent terms. Here we use equation M21

and neglect the integral term equation M22

(The justification for the treatment of these two integrals appears in the Methods section).

Equation (4) exactly matches the equation of radiative transport for a turbid medium with an isotropic phase function, illuminated with a narrow collimated beam, and has been solved by a standard diffusion approximation procedure. This results in the radially dependent diffuse reflectance R_d(ρ) presented in ^{ref. 31},

(6)

where

is the effective attenuation coefficient, equation M25

is the coordinate of the extrapolated boundary condition, equation M26

, and ρ is the exit radius measured from the POE in the surface plane of the semi-infinite medium (see Fig. 1).

Equation (5) describes the phase function correction and is also integrable. The solution yields the phase function correction term, R_p(ρ), for the radially dependent reflectance, obtained by projecting L_p (r,s) on the surface of the medium and integrating over the emerging solid angle

(7)

Physical nature of phase function correction

The physical nature of the phase function correction (7) can be qualitatively understood if we notice thatz + r₁ is the distance a photon would travel through the turbid medium to the detector if it undergoes just a single scattering event, first traveling in the collimated beam to the depth z, then scattering at the angle cos⁻¹ (−z/r₁) and finally traveling distance r₁ to the detector (see Fig. 1). Moreover, the phase function in the integral (−z/r₁) represents the probability of scattering at an angle cos⁻¹. However, note that the attenuation (the exponential in equation (7)) has a characteristic inverse scattering length of equation M28

. The reduced scattering coefficient equation M29

characterizes the attenuation of light due to diffuse scattering, i.e., only diffusely scattered photons are removed from the z + r₁ path; unscattered and low-angle scattered photons remain. If we were considering only unscattered photons, the attenuation factor would have the scattering coefficient μ_s in the exponent, instead of the reduced scattering coefficient equation M30

. The z + r₁ path is not the actual path of real photons, but rather the virtual path of quasi-ballistic photons which are described by equation (7). The contribution of the term describing this behavior can be positive or negative since it depends on the difference of the contributions of the actual phase function and the delta-isotropic phase function. Physically, the PFC term accounts for the contribution of these photons, which have undergone multiple low-angle scattering events plus a single large angle scattering event along the z + r₁ path. That single large angle turn is described by the phase function, thus providing the phase function correction. An illustration of this scattering behavior is shown in Figure 1.

PFC compared with Monte Carlo, δ-P₁ and P₃

To check the accuracy of PFC, we compared the reflectance calculated by Monte Carlo simulation, which is considered to be the gold standard, to the reflectance calculated using PFC, i.e., the sum of the standard DA term, equation (6) and the phase function correction term, equation (7). The Monte Carlo program used was validated against a well accepted code³². Following the procedure described in ^{ref. 33}, we obtained the 95% confidence interval of all the Monte Carlo simulations presented. The margin of error is within 1% of the mean for equation M31

and within 2% for the full range. We include results for two phase functions commonly used in turbid media: the Henyey-Greenstein phase function³⁴ equation M32

and the “Mie” phase function, which can be described using Mie theory³⁵. The Mie phase function was calculated for a Gaussian size distribution with a standard deviation of 10% and a mean size chosen such that the resulting phase function had the desired anisotropy factor, g. The results are presented in Figures 2 and and33.

Figure 2

Sensitivity to the form of the phase function and comparison to other light transport approximations

Figure 3

Effects of absorption and anisotropy factor g

Figure 2a shows the sensitivity of PFC to the particular choice of phase function. While the DA strongly disagrees with the results of the Monte Carlo simulations, especially in the vicinity of the source equation M33

, PFC shows excellent agreement for both the Henyey-Greenstein and Mie phase functions. As an example, a comparison is presented for g 0.95, although for other values of the anisotropy factor g the agreement is equally good. Importantly, both the Henyey-Greenstein and Mie phase functions used in these simulations have identical anisotropy factors yet show marked differences at small source-detector separations. As demonstrated in Figure 2a, these phase function-dependent changes in reflectance are well characterized by PFC. A comparison between PFC and recent attempts to improve the DA is shown in Figure 2b, together with a plot of the percentage error between the approximations and Monte Carlo simulation. The PFC diffusion theory shows dramatically improved agreement with Monte Carlo simulations over both the δ-P₁ and P₃ approximations. For the case of the Mie phase function PFC shows similar results, though the comparison of δ-P₁ and P₃ approximations with Monte Carlo simulations demonstrate improvement compared to Henyey-Greenstein phase function close to the POE (Supplementary Fig. S1).

Figure 3 presents several scenarios where PFC corrects the major deficiencies of the DA. Scenarios of low equation M34

and high

absorption are shown in Figure 3a. The DA gives significant errors for equation M36

, i.e. small distances from the source, however PFC does not suffer from this problem and agrees with the Monte Carlo simulations for all distances. For large source-detector separations and high absorptions, the amount of light emerging from the scattering medium is usually so small that detection is impractical. Figure 3b shows the effects of varying the anisotropy factorg. PFC agrees with Monte Carlo simulations for both g = 0 and g = 0.9 while the DA is only accurate for anisotropy factor g = 0, i.e., totally isotropic scattering. Percentage errors between theory and Monte Carlo for Fig. 3 are shown in Supplementary Fig S2. The inset in Figure 3b depicts the separate contributions of R_p(ρ) and R_d(ρ) and their asymptotic rates of decay for the g = 0.9 case. The contribution of the correction to the reflectance at equation M37

becomes negligible, resulting in an accuracy that is determined by the particular diffusion theory model being used. Note that although we have chosen a well-accepted diffusion theory model to describe isotropic scattering, there may be other solutions to this component of the reflectance.

The correction is intended to work best near the POE where the specific form of the phase function plays a role, while the diffusion approximation is known to work well for large radii where the specific form of the phase function does not matter. The intermediate regime is described less accurately due to the neglected integral term. The average deviation in this range of equation M38

is less than 10%, which is significantly better than the other methods (see Fig. 2b).

Comparison with Monte Carlo simulation within the medium

Though the primary intention of this work is to present a solution to the problem of diffuse reflectance near the POE, we note that the correction also works within the medium. The approach presented above employs the standard diffusion approximation to provide the isotropic solution, which is then corrected with the PFC. However, the standard diffusion approximation is not particularly accurate in predicting the distribution of light near the POE within the medium, even for g = 0. Thus, to check the accuracy of the correction within the medium, we employed Monte Carlo simulation for both isotropic and anisotropic cases using the Henyey-Greenstein phase function and calculated the phase function related change in the fluence rate at several depths within the medium. The PFC fluence rate is obtained by integrating L_p(r,s) over the entire solid angle

(8)

where

. The comparison with Monte Carlo simulations shows good agreement within the medium including the vicinity of the POE (see Fig. 4) at depths of equation M41

. Accurate treatment of larger depths would require including transverse spreading of the collimated beam.

Figure 4

Comparison of the PFC for internal fluence rate with Monte Carlo simulations

Comparison with experimental measurements

We also compared the reflectance calculated using PFC with experimental measurements performed using the setup shown in Fig. 5a (see Methods for details). A high numerical aperture (N.A.) objective was used in this setup because equations (6) and (7) describe a 2π collection solid angle and Monte Carlo simulations show the illumination N.A. has virtually no effect on the shape of the reflectance near the POE (Supplementary Fig. S3). The scattering medium consisted of a mixture of an aqueous glycerin solution and polystyrene microspheres with a diameter of 0.457 μm (Polysciences, Inc.). Using Mie theory³⁵ and assuming a Gaussian size distribution of the microspheres with a standard deviation of 0.011 μm, as listed by the manufacturer, we calculated the reduced scattering coefficient of the scattering sample to be equation M42

mm⁻¹ and the anisotropy factor to be g = 0.89. The absorption coefficient of the scattering sample is approximately μ_a = 0.001 mm⁻¹ (see Supplementary Methods for details). The reflectance signal measured from the polystyrene microsphere sample is presented in Figure 5b. We note that the measurements were performed tens of micrometers from the POE, where the signal strength has significant variation. This distance is comparable to cellular sizes in biological tissue. Thus, measurements of the reflectance from tissue on this scale could be sensitive to alterations in cell architecture associated with normal function or early, pre-cancerous changes.

Figure 5

Experimental schematic and spatially resolved reflectance

Figure 6 shows the experimental reflectance spectra and a comparison with the DA and PFC theories. While the DA exhibits significant deviation from the experiment for distances smaller than 100 μm from the POE, PFC shows excellent agreement.

Figure 6

Comparison with the experiment

Discussion

One of the major applications of this theory will be solving the inverse problem to obtain optical properties of unknown turbid samples especially on a significantly smaller scale than it can be done presently. Fortunately, the same inversion techniques^36,37 that have been used with the standard diffusion approximation to extract equation M43

and μ_a can be applied with the phase function correction, provided that adequate assumptions are made regarding the phase function. Aside from higher accuracy near the source, PFC also allows for the inversion of parameters of the phase function. For example, the anisotropy factor can be measured by assuming either Mie or Henyey-Greenstein scattering and using g as a free parameter in the fitting to experimental data. Moreover, the entire backscattering portion of the phase function can potentially be reconstructed.

A requirement of the technique is the availability of an accurate isotropic solution (g = 0). Though this solution is available for diffuse reflectance, an accurate solution for the fluence rate near the POE is not commonly available. One approach for obtaining the fluence rate within the medium would be to use Monte Carlo simulation for isotropic scattering and combine it with PFC to obtain the fluence for arbitrary phase functions. This is the approach we used for obtaining the results in Figure 6. Thus the result of a single Monte Carlo simulation can be used to obtain solutions for a large array of phase functions and anisotropy parameters. It is also likely that PFC could be used to describe angular dependences in a similar fashion.

In conclusion, PFC, the phase function corrected diffusion theory for light transport in turbid media, addresses the deficiencies of the standard diffusion theory near the POE by accounting for the specific form of the phase function. The technique accurately predicts the correct scattering behavior for two frequently used but very different types of phase functions, the Henyey-Greenstein and Mie theory phase functions (Fig. 2a). The PFC approach results in a substantial improvement over the predictions of the δ-P₁ approximation and the P₃ approximation (Fig. 2b), each of which seeks to improve the performance of the DA. PFC also demonstrates excellent agreement with experimental results (Fig. 6).

We anticipate that this work will be utilized extensively in characterizing photon scattering in turbid media such as human tissue and in a variety of other applications.

Methods

Evaluation of the Integral Terms

To evaluate the first integral term equation M44

we use the standard diffusion approximation expression for L_d (r, s)

(9)

where ϕ_d (r) is the fluence rate and equation M46

is the current is the current density³⁹. By substituting equation M47

in the integral, we get

(10)

From here, we see directly that equation M49

due to the unity normalization of both phase functions.

Let's now simplify the integral term equation M50

and re-write it in a spherical coordinate system (Fig. 7a). The integral can be expressed as μ_s|j_d(r)|∫∫ p (cosθ)(sinθ_j · sin [var phi]

sinθ + cosθ · cosθ_j)sinθdθ_d [var phi]

, where

is the unit vector in the direction of vector j_d(r), θ = c o⁻¹ss ·(s′), and equation M52

The part proportional to the s i nθ_j s i n [var phi]

s · i n θ is zero since equation M53

The part proportional to c o sθ c· o s θ_j gives |j_d (r)| cosθ_j ∫∫ p(cosθ)cosθsinθdθd [var phi]

. Here |j_d(r)|cosθ₀ = s · j_d (r) and, by definition ∫∫p(cosθ)cosθ sin θdθd [var phi]

g. Thus we see that ∫∫p(s, s′)(s′ · j_d (r))ds′ = g (s · j_d (r)). By substituting this into the remaining portion of the integral equation M54

we get

(11)

since both phase functions p(s · s′) and p_DI (s ·s′) have the same anisotropy factor g. Thus indeed μ_s ∫∫ [p(s · s′) − p_DI (s · s′)]L_d (r, s′)ds′ = 0.

Figure 7

Evaluation of the contribution of the neglected integral term

To evaluate the contribution of the second integral term μ_s ∫∫[p(s · s′)−(2π)⁻¹g · δ (1−s · s′)]L_p (r, s′)ds′, we calculated the phase function correction to the radially dependent reflectance including this integral. Then we subtracted R_p (ρ) obtained by neglecting this integral to calculate ΔR_p (ρ). We then plot R_p (ρ) and can see it is just several percent of R_p (ρ for small equation M56

(see Fig. 7b). We also include R_d (ρ) to show that at greater equation M57

both corrections are neglectable. Thus neglecting μ_s ∫∫[p (s · s′)−(2π)⁻¹g · δ(1−s · s′)]L_p (r, s′)ds′ is justified.

Experimental set-up

A schematic of the experimental set-up is shown in Figure 5a. The set-up employs a confocal microscope (FluoView 1000, Olympus) with several modifications. The incident 458 nm laser beam is expanded to fill the back aperture of the objective and passes through a 20/80 beamsplitter (Chroma). The beam is delivered to the scattering sample through a planar apochromatic oil immersion objective (N.A.=1.32, 63× magnification, Olympus). Light scattered from the sample is collected by the same objective and is reflected by a beamsplitter toward the imaging lens and the imaging CCD (DU-434-FI, Andor Technology). The spatial scale for imaging the scattering sample is derived by calibrating the CCD with a counting chamber (Millipore), with a single CCD pixel imaged to 0.327 μm on the sample.

Scattering sample

The sample is placed in an imaging chamber assembled from a glass slide, four adhesive silicone isolators (JTR20R-A2-1.0, Grace Bio-Labs Inc.), and a No. 1 cover glass. The resulting chamber had a 20 mm internal diameter and a 4 mm depth. We placed optical gel (LS-3238, Fiber Optic Center) on the inside surface of the cover glass. After curing, the optical gel formed a layer of approximately 0.1 mm thickness with a refractive index of 1.388. The same refractive index of 1.388 at the 458 nm wavelength is achieved in the scattering sample by mixing a 66% volume fraction of the polystyrene microspheres aqueous solution and a 34% volume fraction of glycerin³⁸. By matching indices we minimize the effects of specular reflection and obtain a matched boundary condition at the interface.

The absorption coefficient of the sample can be estimated as a weighted average of the absorption coefficients of the glycerin, water, and polystyrene in the sample. The absorption coefficient of water⁴⁰ at the wavelength 458 nm is 2.7·10⁻⁵ mm⁻¹. The absorption coefficient of glycerol⁴¹ is approximately 0.003 mm⁻¹ and the absorption coefficient of polystyrene⁴² is less then 1·10⁻⁵ mm⁻¹. Thus the absorption coefficient of the scattering sample is approximately 0.001 mm⁻¹.

Supplementary Material

Supplementary data

Click here to view.^{(91K, doc)}

Acknowledgments

This work was supported by US National Science Foundation grants CBET-0922876 and CBET-0943180, US National Institutes of Health grants R01 EB003472, R01 EB006462 and R33 RR017361, and in part by the Department of Veterans Affairs, Office of Research and Development.

Footnotes

Author contributions E.V. and L.T.P. conceived and developed the theory, L.Q. and L.G. performed the experiments, V.T. and L.Q. performed the data analysis, I.I., E.B.H. and L.T.P. evaluated the method, E.V., V.T., L.Q., I.I., and E.B.H. contributed to the writing of the manuscript, and L.T.P. initiated and supervised the project and wrote the manuscript.

Competing financial interests: The authors declare no competing financial interests.

References

1. Qiu L, et al. Multispectral scanning during endoscopy guides biopsy of dysplasia in Barrett's esophagus. Nature Med. 2010;16:603–606. [PMC free article] [PubMed]

2. Backman V, et al. Detection of preinvasive cancer cells. Nature. 2000;406:35–36. [PubMed]

3. Gurjar RS, et al. Imaging human epithelial properties with polarized light-scattering spectroscopy. Nature Med. 2001;7:1245–1248. [PubMed]

4. Perelman LT, et al. Observation of periodic fine structure in reflectance from biological tissue: A new technique for measuring nuclear size distribution. Phys. Rev. Lett. 1998;80:627–630.

5. Reif R, A'Amar O, Bigio IJ. Analytical model of light reflectance for extraction of the optical properties in small volumes of turbid media. Appl. Opt. 2007;46:7317–7328. [PubMed]

6. Amelink A, Sterenborg H. Measurement of the local optical properties of turbid media by differential path-length spectroscopy. Appl. Opt. 2004;43:3048–3054. [PubMed]

7. Boas DA, Oleary MA, Chance B, Yodh AG. Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media -analytic solution and applications. Proc. Natl Acad. Sci. USA. 1994;91:4887–4891. [PubMed]

8. Itzkan I, et al. Confocal light absorption and scattering spectroscopic microscopy monitors organelles in live cells with no exogenous labels. Proc. Natl Acad. Sci. USA. 2007;104:17255–17260. [PubMed]

9. Chandrasekhar S. Radiative Transfer. Dover; New York: 1960.

10. Yoo KM, Liu F, Alfano RR. When does the diffusion approximation fail to describe photon transport in random media. Phys. Rev. Lett. 1990;64:2647–2650. [PubMed]

11. Blanco S, Fournier R. Short-path statistics and the diffusion approximation. Phys. Rev. Lett. 2006;97:230604-1–230604-4. [PubMed]

12. Chen NG, Bai J. Monte Carlo approach to modeling of boundary conditions for the diffusion equation. Phys. Rev. Lett. 1998;80:5321–5324.

13. Mal'shukov AG, Wang LY, Chu CS, Chao KA. Spin Hall effect on edge magnetization and electric conductance of a 2D semiconductor strip. Phys. Rev. Lett. 2005;95:146601-1–146601-4. [PubMed]

14. vanTiggelen BA, Maynard R, Heiderich A. Anisotropic light diffusion in oriented nematic liquid crystals. Phys. Rev. Lett. 1996;77:639–642. [PubMed]

15. Akkermans E, Wolf PE, Maynard R. Coherent backscattering of light by disordered media - analysis of the peak line-shape. Phys. Rev. Lett. 1986;156:1471–1474. [PubMed]

16. Alerstam E, Svensson T, Andersson-Engels S. Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration. J. Biomed. Opt. 2008;13:060504-1–060504-3. [PubMed]

17. Fang Q, Boas DA. Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units. Opt. Express. 2009;17:20178–20190. [PMC free article] [PubMed]

18. Ren N, et al. GPU-based Monte Carlo simulation for light propagation in complex heterogeneous tissues. Optics Express. 2010;18:6811–6823. [PubMed]

19. Alerstram E, et al. Next-generation acceleration and code optimization for light transport in turbid media using GPUs. Biomed. Opt. Express. 2010;1:658–675. [PMC free article] [PubMed]

20. Yoon G, Prahl SA, Welch AJ. Accuracies of the diffusion approximation and its similarity relations for laser irradiated biological media. Appl. Opt. 1989;28:2250–2255. [PubMed]

21. 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. 1998;E 58:2395–2407.

22. Durian DJ, Rudnick J. Photon migration at short times and distances and in cases of strong absorption. J. Opt. Soc. Am. 1997;A 14:235–245.

23. Kienle A, Patterson MS. Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium. J. Opt. Soc. Am. 1997;A 14:246–254. [PubMed]

24. Seo I, Hayakawa CK, Venugopalan V. Radiative transport in the delta-P-1 approximation for semi-infinite turbid media. Med. Phys. 2008;35:681–693. [PMC free article] [PubMed]

25. Hull EL, Foster TH. Steady-state reflectance spectroscopy in the P-3 approximation. J. Opt. Soc. Am. 2001;A 18:584–599.

26. Menon S, Su Q, Grobe R. Determination of g and mu using multiply scattered light in turbid media. Phys. Rev. Lett. 2005;94:153904-1–153904-4. [PubMed]

27. Chin LCL, Lloyd B, Whelan WM, Vitkin IA. Interstitial point radiance spectroscopy of turbid media. J. Appl. Phys. 2009;105:102025-1–102025-11.

28. Groenhuis RAJ, Ferwerda HA, Tenbosch JJ. Scattering and absorption of turbid materials determined from reflection measurements. 1: Theory. Appl. Opt. 1983;22:2456–2462. [PubMed]

29. Yoon G, Prahl SA, Welch AJ. Accuracies of the diffusion approximation and its similarity relations for laser irradiated biological media. Appl. Opt. 1989;28:2250–2255. [PubMed]

30. Zege EP, Ivanov AP, Katsev IL. Image Transfer Through a Scattering Medium. Springer-Verlag; New York: 1991.

31. Farrell TJ, Patterson MS, Wilson B. A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo. Med. Phys. 1992;19:879–888. [PubMed]

32. Wang LH, Jacques SL, Zheng LQ. MCML - Monte-Carlo modeling of light transport in multilayered tissues. Comput. Methods Programs Biomed. 1995;47:131–146. [PubMed]

33. Liu Q, Ramanujam N. Scaling method for fast Monte Carlo simulation of diffuse reflectance spectra from multilayered turbid media. J. Opt. Soc. Am. 2007;A 24:1011–1025. [PubMed]

34. Henyey LG, Greenstein JL. Diffuse radiation in the galaxy. Astrophys. J. 1941;93:70–83.

35. Toublanc D. Henyey-Greenstein and Mie phase functions in Monte Carlo radiative transfer computations. Appl. Opt. 1996;35:3270–3274. [PubMed]

36. Zonios G, et al. Diffuse reflectance spectroscopy of human adenomatous colon polyps in vivo. Appl. Opt. 1999;38:6628–6637. [PubMed]

37. Sun J, et al. Influence of fiber optic probe geometry on the applicability of inverse models of tissue reflectance spectroscopy: computational models and experimental measurements. Appl. Opt. 2006;45:8152–8162. [PubMed]

38. Fang H, et al. Confocal light absorption and scattering spectroscopic microscopy. Appl. Opt. 2007;46:1760–1769. [PubMed]

39. Wang LV, Wu HI. Biomedical Optics: Principles and Imaging. Wiley; New Jersey: 2007.

40. Hale GM, Querry MR. Optical constants of water in the 200-nm to 200-microm wavelength region. Appl. Opt. 1973;12:555–563. [PubMed]

41. Reich O, Lohmannsroben H-G, Schael F. Optical sensing with photon density waves: Investigation of model media. Phys. Chem. Chem. Phys. 2003;5:5182–5187.

42. Firbank M, Arridge SR, Schweiger M, Delpy DT. An investigation of light transport through scattering bodies with non-scattering regions. Phys. Med. Biol. 1996;41:767–783. [PubMed]