Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Nat Commun. Author manuscript; available in PMC 2012 December 13.
Published in final edited form as:
PMCID: PMC3370954

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


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 tissue18. 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 generations912. Propagation of electrons in semiconductors13, anisotropic light diffusion in liquid crystals14, and coherent backscattering15 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,

s[center dot][nabla]L(r,s)=(μa+μs)L(r,s)+μsp(s,s)L(r,s)ds+S(r,s)

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)9.

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 time1619 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 isotropic20. It is not sensitive to the particular form of the phase function, reducing phase function information to a single parameter, g=p(s[center dot]s)(s[center dot]s)ds, called the average cosine of the phase function or the anisotropy factor, which in turn is incorporated into the reduced scattering coefficient, μs=μs(1g). Hence, μs 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 functions2123. One such approach, called the δ-P1 approximation, is to solve the equation for the delta-Eddington phase function, which results in some improvements for characterizing forward scattering media24. 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 P3 approximation uses only the 0th through 3rd terms in the Legendre series25. 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 used26,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 1/μs) 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.


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

p(s[center dot]s)=pDI(s[center dot]s)+[p(s[center dot]s)pDI(s[center dot]s)].

where pDI(s · s′)=[1− g +2g ·†(1−s·s′)]/4π is the delta-isotropic phase function28. 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

S[center dot][nabla]L(r,s)=(μa+μs)L(r,s)+μs4πL(r,s)ds+μs[p(s[center dot]s)PDI(s[center dot]s)]L(r,s)ds.

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, Ld(r, s) a phase function dependent part, (Lp(r, s) and an on-axis part due to the narrow collimated beam, Lc(r, s).

After substitution for L(r, s) in equation (3), we seek an expression for Lc(r, s) which satisfies the boundary condition for a narrow collimated beam and a portion of equation (3), s[center dot][nabla]Lc(r,s)=(μa+μs)Lc(r,s). Note that this portion is described by μs rather than μs and specifically accounts for low-angle scattered photons29 (see Fig. 1 for illustration). The boundary condition corresponding to a narrow collimated beam entering a semi-infinite medium is L[mid ]z=0=P0(2π)2ρ1δ(ρ)δ(1s[center dot]z^) for inward directions (s[center dot]z^>0), where p0 is the power of the incoming beam, ρ and z are the cylindrical coordinates and z^ is a unit vector along the z axis. This yields30 Lc(r,s)=P0(2π)2ρ1δ(ρ)δ(1s[center dot]z^)exp[(μa+μs)z] and the boundary conditions Ld[mid ]z=0=0 and Lp[mid ]z=0 for s[center dot]z^>0.

The transport equation (3) can then be written as a system of two equations for the remaining parts of the radiance Ld(r, s) and Lp(r, s). It is natural then to distribute the two remaining Lc(r, s) dependent terms so that the isotropic, phase function independent term μs4πLc(r,s)ds is combined is combined with the equation for the predominately isotropic Ld(r, s) and the anisotropic, phase function dependent term μs[p(s[center dot]s)pDI(s[center dot]s)]Lc(r,s)ds is combined with the equation for the anisotropic Lp(r, s).

The first equation,

s[center dot][nabla]Ld(r,s)=(μa+μs)Ld(r,s)+μs4πLd(r,s)ds+P0μs8π2δ(ρ)ρexp[(μa+μs)z]

leads to the standard diffusion approximation and the second equation,

s[center dot][nabla]Lp(r,s)=(μs+μs)Lp(r,s)+P0μs2πδ(ρ)ρ[p(s[center dot]z^)1g4π]exp[(μs+μs)z],

retains all of the phase function dependent terms. Here we use μs[p(s[center dot]s)pDI(s[center dot]s)]Ld(r,s)ds=0 and neglect the integral term μs[p(s[center dot]s)(2π)1g[center dot]δ(1s[center dot]s)]Lp(r,s)ds (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 Rd(ρ) presented in ref. 31,


where μeff=[3μa(μa+μs)]1/2 is the effective attenuation coefficient, zb=2/[3(μa+μs)] is the coordinate of the extrapolated boundary condition, r1=ρ2+z2,r2=ρ2+(z+2zb)2, 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, Rp(ρ), for the radially dependent reflectance, obtained by projecting Lp (r,s) on the surface of the medium and integrating over the emerging solid angle



Physical nature of phase function correction

The physical nature of the phase function correction (7) can be qualitatively understood if we notice thatz + r1 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−1 (−z/r1) and finally traveling distance r1 to the detector (see Fig. 1). Moreover, the phase function in the integral (−z/r1) represents the probability of scattering at an angle cos−1. However, note that the attenuation (the exponential in equation (7)) has a characteristic inverse scattering length of μs. The reduced scattering coefficient μs characterizes the attenuation of light due to diffuse scattering, i.e., only diffusely scattered photons are removed from the z + r1 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 μs. The z + r1 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 + r1 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, δ-P1 and P3

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 code32. 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 ρμs<0.1 and within 2% for the full range. We include results for two phase functions commonly used in turbid media: the Henyey-Greenstein phase function34 pHG(s[center dot]s)=[(1g2)/4π]/[1+g22g(s[center dot]s2)]3/2 and the “Mie” phase function, which can be described using Mie theory35. 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 (ρμs<0.1), 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 δ-P1 and P3 approximations. For the case of the Mie phase function PFC shows similar results, though the comparison of δ-P1 and P3 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 (μa/μs=0.01) and high (μa/μs=1) absorption are shown in Figure 3a. The DA gives significant errors for ρμs<0.1, 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 Rp(ρ) and Rd(ρ) and their asymptotic rates of decay for the g = 0.9 case. The contribution of the correction to the reflectance at ρμs>>1 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 ρμs[open square]0.21 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 Lp(r,s) over the entire solid angle


where ρ1=ρ2+(zz)2. The comparison with Monte Carlo simulations shows good agreement within the medium including the vicinity of the POE (see Fig. 4) at depths of zμs0.5. 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 theory35 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 μs=5.6 mm−1 and the anisotropy factor to be g = 0.89. The absorption coefficient of the scattering sample is approximately μa = 0.001 mm−1 (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


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 techniques36,37 that have been used with the standard diffusion approximation to extract μs 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 δ-P1 approximation and the P3 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.


Evaluation of the Integral Terms

To evaluate the first integral term μs[p(s[center dot]s)pDI(s[center dot]s)]Ld(r,s)ds we use the standard diffusion approximation expression for Ld (r, s)

Ld(r,s)=14π[ϕd(r)3s[center dot]jd(r)]

where ϕd (r) is the fluence rate and jd(r)=13(μa+μs)[nabla]ϕd(r) is the current is the current density39. By substituting 14π[ϕd(r)3s[center dot]jd(r)] in the integral, we get

μs[p(s[center dot]s)pDI(s[center dot]s)]Ld(r,s)ds==μs4π[p(s[center dot]s)pDI(s[center dot]s)]ϕd(r)ds3μs4π[p(s[center dot]s)pDI(s[center dot]s)](sjd(r))ds

From here, we see directly that μs4π[p(s[center dot]s)pDI(s[center dot]s)]ϕd(r)ds=μsϕd(r)4π[p(s[center dot]s)pDI(s[center dot]s)]ds=0 due to the unity normalization of both phase functions.

Let's now simplify the integral term p(s,s)(s[center dot]jd(r))ds and re-write it in a spherical coordinate system (Fig. 7a). The integral can be expressed as μs|jd(r)|∫∫ p (cosθ)(sinθj · sin[var phi] sinθ + cosθ · cosθj)sinθdθd[var phi], where j^ is the unit vector in the direction of vector jd(r), θ = c o−1ss ·(s′), and θj=co1sj^d([center dot]s) The part proportional to the s i nθj s i n[var phi] s · i n θ is zero since 02πs[var phi].d[var phi]=0 The part proportional to c o sθ c· o s θj gives |jd (r)| cosθj ∫∫ p(cosθ)cosθsinθdθd[var phi]. Here |jd(r)|cosθ0 = s · jd (r) and, by definition ∫∫p(cosθ)cosθ sin θdθd[var phi]g. Thus we see that ∫∫p(s, s′)(s′ · jd (r))ds′ = g (s · jd (r)). By substituting this into the remaining portion of the integral 3μs4π[p(s[center dot]s)pDI(s[center dot]s)](s[center dot]jd(r))ds we get

3μs4π[p(s[center dot]s)pDI(s[center dot]s)](s[center dot]jd(r))ds=3μs4π[g(s[center dot]jd(r))g(s[center dot]jd(r))]=0,

since both phase functions p(s · s′) and pDI (s ·s′) have the same anisotropy factor g. Thus indeed μs ∫∫ [p(s · s′) − pDI (s · s′)]Ld (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π)−1 g · δ (1−s · s′)]Lp (r, s′)ds′, we calculated the phase function correction to the radially dependent reflectance including this integral. Then we subtracted Rp (ρ) obtained by neglecting this integral to calculate ΔRp (ρ). We then plot Rp (ρ) and can see it is just several percent of Rp (ρ for small ρμs (see Fig. 7b). We also include Rd (ρ) to show that at greater ρμs both corrections are neglectable. Thus neglecting μs ∫∫[p (s · s′)−(2π)−1 g · δ(1−s · s′)]Lp (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 glycerin38. 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 water40 at the wavelength 458 nm is 2.7·10−5 mm−1. The absorption coefficient of glycerol41 is approximately 0.003 mm−1 and the absorption coefficient of polystyrene42 is less then 1·10−5 mm−1. Thus the absorption coefficient of the scattering sample is approximately 0.001 mm−1.

Supplementary Material

Supplementary data


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.


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.


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]