PROCEEDINGS OF SPIE
SPIEDigitalLibrary.org/conference-proceedings-of-spie
Diffusion theory for multilayered
scattering media
Joseph Hollmann, Lihong V. Wang
Joseph Hollmann, Lihong V. Wang, "Diffusion theory for multilayered
scattering media," Proc. SPIE 5695, Optical Interactions with Tissue and Cells
XVI, (15 April 2005); doi: 10.1117/12.590886
Event: SPIE BiOS, 2005, San Jose, CA, United States
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/25/2018 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
Diffusion theory for multi-layered scattering media
Joseph Hollmann and Lihong V. Wang
Optical Imaging Laboratory, Department of Biomedical Engineering, Texas A&M
University
3120 TAMU, College Station, Texas 77843-3120, USA
ABSTRACT
We present an approximation to light propagation through skin using an improved multilayered diffusion
equation. This equation allows for the influence of multiple layers on diffusely reflected light to be
considered; which is accomplished by embedding imaginary point sources in the medium that correspond
to each layer’s transport mean free path. The new approximation is then compared to a previous two-layer
solution by Farrell et al. to verify improvement as well as to results from Monte Carlo simulations to verify
accuracy.
Keywords
: tissue,
Optical diagnostics for medicine, Medical optics and biotechnology, Spectroscopy,
Diffusion
1. INTRODUCTION
There is great interest in determining the optical parameters of biological tissue non-invasively for a wide
range of applications in therapeutics and diagnostics. The first step in this process is to develop an accurate
and relatively fast method that can be used to simulate the propagation of light through tissue. This is
hampered by the difficulty that light is highly scattered by most tissue. This means that the optical signal
undergoes a large number of scattering events before its photons exit the medium. A signal that contains
photons with no correlation between their entering and exiting directions are likely to have gone through
large distances in the tissue. While these photons carry information about the medium, their chaotic paths
are difficult to resolve. The challenge is in developing a forward method to approximate the signal inside
the tissue.
One technique for finding the diffuse reflectance from highly scattering tissue is to use the homogenous
diffusion approximation to the Radiative Transport Equation. This method assumes that the light reflected
from tissue due to an incident collimated beam can be modeled as light coming from an equivalent
isotropic point source located within a homogenous medium. For this to work, the tissue’s reduced
scattering (
μ
s
’) is assumed to be much higher than the absorption coefficient (
μ
a
). The reduced scattering
coefficient is arrived at by
)
1
(
'
g
s
s
−
=
μ
μ
(1)
where
μ
s
and g are the scattering and anisotropy coefficients respectively. The isotropic point source is
located at one transport mean free path inside the medium (
/
1
μ
t
’), where
a
s
t
μ
μ
μ
+
=
'
'
(2)
Detailed discussions of this method for normal incidence can be found in Farrell et al.
1
and Oblique
incidence in Wang et al.
2
As mentioned above these two methods model tissue as a homogenous medium. For skin and other
complicated tissues this is not an accurate assumption. At a minimum skin contains both the epidermis and
dermis layers. If these are optically thin we also have to take into account how the subcutaneous tissue
influences the light signal. To improve upon this assumption Kienle et al.
3
have developed a two-layer
Optical Interactions with Tissue and Cells XVI, edited by Steven L. Jacques,
William P. Roach, Proc. of SPIE Vol. 5695 (SPIE, Bellingham, WA, 2005)
1605-7422/05/$15 · doi: 10.1117/12.590886
101
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/25/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
solution to the diffusion equation that models the light reflected from skin as coming from a single
equivalent point source. The solution assumes that the first layer is optically thick and the isotropic point
source is located at one transport mean free path within it. This solution does not fully incorporate the
effects of the second layer; which is a problem if the first is optically thin. We present here a solution to
the diffusion equation that utilizes an equivalent point source for each of the two layers; which allows us to
overcome these difficulties.
2. METHODOLOGY
2.1. Theory
The problem of interest has an infinitesimal, collimated beam incident at a normal angle onto scattering
medium made up of two homogenous layers. The coordinate system is set up so that its origin is located at
the point where the beam enters the medium with its z+ axis is pointing down. The first layer has a
thickness of L and the second layer extends to infinity. Both layers are considered to be radially
homogenous. The resulting fluence rate (
Φ
) and reflectance is modeled by replacing the incident
collimated beam with an equivalent isotropic point source in each layer. The geometry of this problem is
shown in Figure 1. The optical properties for layers one and two are denoted by their subscripts.
The locations of the point sources (z
0i
for i=1,2) are given by the mean intensity over the total intensity of
light’s fluence within each layer.
Figure 1: Geometry of two-
layer solution to the diffusion equation. The positions of the equivalent isotropic point
sources are denoted by z
01
and z
02
. L is the thickness of layer 1 and the positive z direction is into the tissue.
102 Proc. of SPIE Vol. 5695
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/25/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
'
2
'
2
'
1
'
2
'
1
02
'
1
'
1
'
1
'
1
'
1
0
'
1
0
01
1
))
(
(
exp
))
(
(
exp
))
exp(
1
(
)
exp(
)
1
(
1
)
(
exp
)
(
exp
t
t
t
L
t
t
L
t
t
t
t
t
L
t
L
L
dz
L
z
L
dz
L
z
L
z
z
L
L
L
dz
L
dz
L
z
z
μ
μ
μ
μ
μ
μ
μ
μ
μ
μ
μ
+
=
−
−
−
−
−
−
=
−
−
−
+
−
=
−
−
=
∫
∫
∫
∫
∞
∞
(3)
where L is the thickness of the first layer. The corresponding source strengths (or weightings) for these
positions are
∫
∫
∞
−
=
−
−
−
=
−
−
=
−
=
L
t
t
t
t
w
L
t
t
t
w
L
a
dz
L
L
z
a
z
S
L
a
dz
z
a
z
S
)
exp(
)
)
(
exp(
)
(
))
exp(
1
(
)
exp(
)
(
'
1
'
2
'
1
'
2
'
2
'
2
2
0
'
1
'
1
'
1
'
1
'
1
1
μ
μ
μ
μ
μ
μ
μ
(4)
where
a
i
’
is the transport albedo for each layer:
a
i
’=
μ
si
’/
μ
ti
’,
i=1 and 2. We can see that if the first layer
becomes too thick,
S
w2
approaches zeros and z
01
approaches 1/
μ
t1
’. This then becomes the single-layer
diffusion equation.
The steady-state diffusion equation becomes
2
11
11
1
01
2
2
222
2
02
(
)() (,,)
(
)() (,,
)
aw
aw
Dr
rSxyzz
D
rrSxyzz
μδ
μδ
∇Φ
− Φ
=−
−
∇Φ
−
Φ
=−
−
(5)
where
)
(
3
1
'
si
ai
i
D
μ
μ
+
=
is the diffusion constant and
r=(x,y,z)
. The dirac delta functions correspond to
each source position.
These equations can be solved for by first taking their two-dimensional Fourier transform:
)
(
)
,
,
(
)
,
,
(
0
2
2
2
i
i
wi
b
a
i
i
b
a
i
z
z
D
S
s
s
z
s
s
z
z
−
−
=
−
∂
∂
δ
φ
α
φ
i
=1,2
(6)
where
s
a
and
s
b
are spatial frequency dummy variables. We can take advantage of the radial homogeneity
by saying that
s
i
2
=
s
ai
2
+
s
bi
2
and by introducing the variable
i
ai
i
i
i
D
s
D
μ
α
+
=
2
2
.
Using a Green’s function and the following boundary conditions can solve for equation 4:
First, we use the extrapolated boundary condition for a matched index of refraction at the tissue-air
interface to give us
Proc. of SPIE Vol. 5695 103
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/25/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
0
)
,
(
1
=
−
s
z
b
φ
(7)
where
1
2
b
z
AD
=
(8)
is the position of the extrapolated boundary. The negative sign implies the boundary is located a
bove the
medium. The boundary condition is handled in the term A
1
1
e
ff
e
ff
r
A
r
+
=
−
(9)
where r
eff
is the number of photons internally, diffusely reflected at the boundary. For an index-matched
boundary, r
eff
is zero and A=1. Haskell et al.
4
discusses how to find r
eff
for index mismatches at the
boundary.
Using the supposition th
at the second layer is infinite we get
0
)
,
(
2
=
∞
s
φ
(10)
If we assume tissue layers have similar indices of refraction, we get
1
)
,
(
)
,
(
2
2
2
1
2
1
=
=
n
n
s
L
s
L
φ
φ
(11)
L
z
L
z
z
s
z
D
z
s
z
D
=
=
∂
∂
=
∂
∂
)
,
(
)
,
(
2
2
1
1
φ
φ
(12)
We then get the solution to equation 4 to be
()
22
1
1
1
01
1
01
2
2
02
11
1
22
1
11
1
22
11
1012
202
1
11
exp(
(
))
cosh( (
))
sinh( (
))
exp(
(
))
(,)
2
sinh( (
))
cosh( (
))
exp(
)
1
sinh( (
))
exp(
(
)
bw
w
bb
wbw
D
z z S
Lz
Lz
S
z L
D
zs
DzLDzL
D
LS
z z
S
z
L
z
D
α
αααα
α
φ
αα
αα
α
αααα
α
−+
−−
− −
− −
=
++
+
−−
++−−+
+
()
1
1101
22
1
11
1
11
) exp(
)
exp(
)
2
sinh( (
))
cosh( (
))
2
b
w
bb
z
S
zz
DzLDzL
D
α
α
αα αα
α
−−
+
++
+
(13)
for 0<z<L and
104 Proc. of SPIE Vol. 5695
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/25/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
()
()
()
()()
()
11012
2
222
1
11
1
11
1
2
2
2
02
1
1
01
22
222
1
11
1
1
exp
(
)
(, )
2 exp(
)
sinh( (
))
cosh( (
))
c
osh
(
) exp(
)
exp
(
)
exp
(
)
2 exp(
)
sinh( (
))
cosh( (
))
sinh
wb
bb
bww
bb
Szzz
zs
LD
zL D
zL
D
Lz
z
S
z L S
Lz
D
LD
zL D
zL
αα
φ
αα
α
α
α
α
αα
α
α
α
αα
α
α
α
α
−−+−
=
−+++
+−
−−−−−
−
−+++
−
()
( )( )
()
()
()
22
202
1
101
222
1
11
1
2201
22
(
) exp(
)
exp
(
)
exp
(
)
2 exp(
)
sinh( (
))
cosh( (
))
exp
2
bw
w
bb
w
Lz
zS
z L S
Lz
LD
zL D
zL
Szz
D
αα
α
αα
α
α
α
α
α
+−
−−+−−
−+++
−−
+
(14)
for z > L.
We then need to take the inverse Fourier transform to convert equations 13 and 14 into their final solutions.
∫∫
∞
∞
−
∞
∞
−
+
−
=
Φ
2
1
2
1
2
))
(
exp(
)
,
(
4
1
)
,
,
(
ds
ds
y
s
x
s
i
s
z
z
y
x
i
i
φ
π
(15)
To simplify equation 15 we can take advantage of the radial homogeneity of the model and convert to
cylindrical coordinates.
ds
s
J
s
z
s
z
i
i
∫
∞
=
Φ
0
0
)
(
)
,
(
2
1
)
,
(
ρ
φ
π
ρ
(16)
where J
0
is a zeroth order Bessel function of the first kind. These equations cannot be solved in closed
form, so we employ the Simpson’s quadrature
5
to find the numerical solution. It is not necessary to
integrate to infinity and infact, we must impose a limiting condition due to the hyperbolic functions or else
computer round off error could drastically influence the final answer. For this situation, it can be assumed
that the function
φ
i
is monotonically decreasing as z grows. If
φ
i
changes direction and begins to grow, we
know that the error is influencing the results and should end our integration.
For an index-matched boundary, the spatially resolved reflectance is given from the fluence as
6
1
11
0
(
,)
11
()
(, 0)
42
z
z
RzD
z
ρ
ρρ
=
∂Φ
=Φ
= +
∂
(17)
this equation can be simplified by using
4
1
1
0
(, 0)
z
zA
z
ρ
=
∂Φ
Φ==
∂
(18)
Substituting equation 18 into 17 and using A=1 from a
bove, we get
Proc. of SPIE Vol. 5695 105
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/25/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
1
1
(
)(,
0)
2
Rz
ρρ
=Φ
=
(19)
This discussion has been limited to a two-layer tissue structure but can be extended to any number of
layers. First place an equivalent isotropic point source in each layer using equations 3 and 4 for the
locations and weightings. Then use equations 11 and 12 to find the boundary conditions at the borders of
each layer and then solve for the constants of integration.
2.2. Simulations
For the purposes of this paper the Monte Carlo simulation shall be the “gold standard.” A detailed
discussion can be found in Wang et all
7
but for our purposes it will suffice to say that Monte Carlo uses a
weighted random walk to predict the paths of individual photons through a turbid medium. Since the
simulation is statistical in nature, it requires a large number of photons to be run so that the results are
within acceptable error bounds. The two-source solution presented within this paper and the single-source
solution mentioned above
3
were then programmed so that we could verify the improved accuracy of our
method.
Two Monte Carlo simulations were run to test the approximations using the parameters listed in Figure 2.
The optical parameters were kept the same for each simulation with only the depth changing. This is to test
the validity of each approximation for an optically thick and thin first layer. Optically thick is defined as
L>1/
μ
t1
’.
Run 1
Run 2
μ
s1
’ (cm-1)
10 10
μ
s2
(cm-1)
12 12
μ
a1
(cm-1)
.2 .2
μ
a2
(cm-1)
.01 .01
L
(cm)
.12 .05
Figure 2: Optical Parameters for simulations
3. RESULTS
Figure 3: Spatially resolved reflectance curves for the Monte Carlo Simulation (symbols), and the
single (solid)
and double source (dotted) approximations.
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
10
0
10
1
reflectance [W/cm
2
]
distance [cm]
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
10
0
10
1
reflectance [W/cm
2
]
distance [cm]
106 Proc. of SPIE Vol. 5695
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/25/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
The results for the reflectance due to the parameters listed under run 1 are shown in Figure 3. The
reflectance from the two-source approximation seems to agree with the Monte Carlo results at a distance of
0.1 cm from the source. We can see from figure 4 that the two-source approximation’s error does not
exceed 7%. We can see that the single source approximation has a maximum error of about 10% in the
same area.
We can see that the reflectance due to the two-source diffusion approximation becomes accurate much
earlier - at a distance of less than 0.05 cm. It is not generally true that the two-source approximation is
accurate closer to the source than the single source approximation.
The results for an optically thin layer (Figure 2, run 2) are shown in Figure 5. We see that the reflectance
calculated from the single-source approximation oscillates wildly as compared to the Monte Carlo results.
The error is above
100% at the local maxima’s and minima’s of the single-source reflectance curve.
Figure 5: Spatially resolved reflectance curves for the Monte Carlo Simulation (symbols), and the single
(solid) and double source (dotted) approximations
Figure 4: Error of the single (circles) and double-
source (squares) approximations
versus Monte Carlo results
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
-0.15
-0.1
-0.05
0
0.05
0.1
error [-]
distance [cm]
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
-0.15
-0.1
-0.05
0
0.05
0.1
error [-]
distance [cm]
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
10
-1
10
0
10
1
10
2
reflectance [W/cm
2
]
distance [cm]
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
10
-1
10
0
10
1
10
2
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
10
-1
10
0
10
1
10
2
reflectance [W/cm
2
]
distance [cm]
Proc. of SPIE Vol. 5695 107
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/25/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
Clearly the approximation fails for in this instance. This is characteristic of the single-source solution when
the first layer’s thickness is less than one transport mean free path. This is because the single-source
approximation does not take into account the second layer’s contribution to the reflected light.
The two-source approximation reflectance closely follows the Monte Carlo results for an optically thin
layer with and in figure 6, we see that it has an error of less than 5%. even though the previous
approximation fails. This means that the two-source solution effectively reduces our error’s dependency on
the layer depth.
4. CONCLUSIONS
This paper has laid out an improved method for solving the diffusion approximation for light traveling
through multiple layers. The method was illustrated by deriving the diffusion approximation for the
specific case of a two-layer, turbid medium. A collimated, incident beam was replaced by one equivalent
isotropic point in each layer. These sources’ locations and strengths were based on their layer’s optical
properties. The resulting diffusion equations were then solved by using a two-dimensional Fourier
transform to get them in to ordinary differential equation form. We then used the Simpson’s quadrature
take the inverse Fourier transform and get the fluence of light through the medium. Then using equation
15, we converted the fluence at z=0 into reflected light. The results were then compared against the results
of previous solutions. We also briefly addressed methods for extending the solution to x layers using x
sources, where x is a positive, finite integer.
Figure 3 shows both the approximations agree with the Monte Carlo simulations for an optically thick first
layer. We also see that the accuracy region of the two-source approximation is dependent on both layers
where as the single-source’s solution only depends on the top layer. In Figure 5 we see that the previous
solution completely breaks down for layer thicknesses less than 1/
μ
t1
’ while the two-source approximation
stays within a 5% error bound.
Figure 6: Error of the single (circles) and double-
source (squares) approximations versus Monte
Carlo results
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
-1.5
-1
-0.5
0
0.5
1
1.5
reflectance [W/cm
2
]
distance [cm]
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
-1.5
-1
-0.5
0
0.5
1
1.5
reflectance [W/cm
2
]
distance [cm]
108 Proc. of SPIE Vol. 5695
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/25/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
REFERENCES
1. T. .J. Farrell and M. S. Patterson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance
for the noninvasive determination of tissue optical properties
in vivo
,” Med. Phys. 19, 879 - 888 (1972)
2.
L.-H. Wang and S. L. Jacques, "Use of a laser beam with an oblique angle of incidence to measure the
reduced scattering coefficient of a turbid medium," Applied Optics 34, 2362-2366 (1995).
3. A. Kienle, M.S.Patterson, N.Dognitz, R.Bays, G.Wagnieres, and H. Van Den Bergh, “Noninvasive determination of
the optical properties of two-layered turbid media,” Appl. Opt 37, 779-791 (1998).
4. R.C. Haskell, L.O. Svaasand, T.T. Tsay, T.C. Feng, M. McAdams, and B.J. Tromberg, ”Boundary Cond
itions for
the diffusion equation in Radiative Transfer,” J. Opt. Soc. Am A, 11, 3737-3741 (1994).
5. M. N. O. Sadiku, “Finite Difference Methods,” in
Numerical Techniques in Electromagnetics, 2
nd
edition
(CRC
Press, Boca Raton, Florida, 2001) pp 197-199.
6
6. G. Alexandrakis, T.J. Farrell, M.S. Patterson, ”Accuracy of the diffusion approximation in determining
the optical properties of a two-layer turbid medium,” Appl. Opt. 37, 7401-7409 (1998)
7. L.-H. Wang, S. L. Jacques, and L.-Q. Zheng, "MCML - Monte Carlo modeling of photon transport in multi-
layered tissues," Computer Methods and Programs in Biomedicine 47, 131-146 (1995).
Proc. of SPIE Vol. 5695 109
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/25/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use