628
J. Opt. Soc. Am. A / Vol. 31, No. 3 / March 2014
C. Sandoval and A. D. Kim
Deriving Kubelka–Munk theory from radiative transport Christopher Sandoval and Arnold D. Kim* Applied Mathematics Unit, School of Natural Sciences, University of California, Merced, 5200 North Lake Road, Merced, California 95343, USA *Corresponding author: [email protected] Received November 12, 2013; accepted January 6, 2014; posted January 16, 2014 (Doc. ID 201164); published February 21, 2014 We derive Kubelka–Munk (KM) theory systematically from the radiative transport equation (RTE) by analyzing the system of equations resulting from applying the double spherical harmonics method of order one and transforming that system into one governing the positive- and negative-going fluxes. Through this derivation, we establish the theoretical basis of KM theory, identify all parameters, and determine its range of validity. Moreover, we are able to generalize KM theory to take into account general boundary sources and nonhomogeneous terms, for example. The generalized Kubelka–Munk (gKM) equations are also much more accurate at approximating the solution of the RTE. We validate this theory through comparison with numerical solutions of the RTE. © 2014 Optical Society of America OCIS codes: (290.4210) Multiple scattering; (290.7050) Turbid media; (030.5620) Radiative transfer; (000.3860) Mathematical methods in physics. http://dx.doi.org/10.1364/JOSAA.31.000628
1. INTRODUCTION Multiple scattering of light in a turbid medium is governed by the theory of radiative transport [1,2]. This theory provides a complete description of absorption, scattering, and radiation of light. However, applying radiative transport theory is challenging for practical problems because exact solutions of the radiative transport equation (RTE) are known only for simple problems. Even numerical solutions are challenging due to the large number of variables involved. Kubelka and Munk [3,4] developed a phenomenological theory that describes the flow of power in positive- and negative-going directions, which we denote by F , respectively, through a plane parallel slab of a scattering medium. The Kubelka–Munk (KM) equations are dF K SF SF − ; dτ
(1.1a)
dF − K SF − SF ; dτ
(1.1b)
−
where τ is the optical depth, and K and S are related to absorption and scattering, respectively. This so-called two-flux theory is intuitive and much simpler than the RTE. Hence, it has been applied to several practical problems (see [5] and references contained therein). However, the theoretical basis of KM theory is not entirely established [2]. Several studies [6–18] have addressed the theoretical basis for KM theory and its applicability to interpret measured data. However, a systematic derivation of the KM equations from the RTE is still needed to determine all coefficients, to determine its range of validity, and to provide a framework to derive extensions that will allow for broader application to practical problems. We present here a systematic derivation of the KM equations from the RTE. This derivation comes from analyzing 1084-7529/14/030628-09$15.00/0
the system of equations resulting from applying the double spherical harmonics method of order one to a general boundary value problem for the one dimensional RTE. This theoretical framework also provides a means to generalize KM theory to a larger class of problems including general boundary conditions and nonhomogeneous terms. For example, we identify the fundamental reason why a four-flux, rather than a two-flux, theory is needed to take into account collimated sources. We validate all aspects of this theory through comparisons with numerical solutions of the RTE. A key step in this derivation is the application of the double spherical harmonics method [19–22] to solve the general boundary value problem for the one dimensional RTE. The double spherical harmonics method uses two different expansions in the angular variable to represent the specific intensity over two distinct hemispheres. Doing so also provides a natural representation for the solution which, in general, is discontinuous at directions tangent to the parallel boundary planes. The double spherical harmonics method provides a natural way to prescribe boundary conditions since they are prescribed over half ranges. It is this piecewise representation of the specific intensity that allows for the systematic derivation of the KM equations governing the positive- and negativegoing fluxes. The remainder of this paper is as follows. In Section 2, we present the general boundary value problem for the one dimensional RTE. In Section 3, we apply the double spherical harmonics method to solve this boundary value problem for the RTE. In Section 4, we transform the double spherical harmonics system of order one to a 4 × 4 system, governing the positive- and negative-going fluxes and two auxiliary functions, which we call the generalized Kubelka–Munk (gKM) equations. In Section 5, we reduce this generalized system to the 2 × 2 system of KM equations, given in Eq. (1.1), assuming that multiple scattering is strong. We give numerical © 2014 Optical Society of America
C. Sandoval and A. D. Kim
Vol. 31, No. 3 / March 2014 / J. Opt. Soc. Am. A
results of comparisons between results of this theory and those of direct numerical solution of the RTE in Section 6. In Section 7, we give our conclusions.
−μ
∂I − ϖ I− − 0 ∂τ 2
1
−1
hμ; μ0 Iμ0 ; τdμ0 Qμ; τ in τ1 < τ < τ2 ; (2.1a)
Iμ; τ1 f μ on 0 < μ ≤ 1;
(2.1b)
Iμ; τ2 f − μ on − 1 ≤ μ < 0:
(2.1c)
Here, Iμ; τ is the specific intensity which depends on μ cos θ, where θ is the polar angle, and τ is the optical depth. We denote the albedo by ϖ 0 , which is defined in terms of the absorption and scattering coefficients, μa and μs , respectively, as ϖ 0 μs ∕μs μa . The function Qμ; τ is a source, and the functions f give sources on the boundaries, τ τ1 and τ τ2 , respectively, for all directions pointing into the plane-parallel slab. The redistribution function, h, gives the fraction of light scattered in direction μ due to light incident in direction μ0 . It is normalized according to Z
1 2
1
−1
hμ; μ0 dμ0 1:
(2.2)
Z
1
−1
hμ; μ0 μ0 dμ0 gμ;
(2.3)
(3.1c)
I − μ; τ2 f − −μ;
(3.1d)
Here, we have assumed that scattering is spherically symmetric, so that h−μ; −μ0 hμ; μ0 and hμ; −μ0 h−μ; μ0 . We p introduce the functions P~ n μ 2n 1P n 2μ − 1 defined over the half-range, 0 < μ ≤ 1, where fP n g is the set of Legendre polynomials. We call fP~ n g the set of normalized, half-range Legendre polynomials. They satisfy the following orthogonality relation: Z
1
I μ; τ
∞ X
~ c n τP n μ:
(3.3)
n0
Let c denote the vectors whose components are c T c 0 ; c1 ; . By substituting Eq. (3.3) and making use of orthogonality relation (3.2), we find that c satisfies the system d M dτ 0
0 −M
c c−
ϖ 0 H 1 − I ϖ 0 H 2 q ; q−
ϖ 0 H 2 ϖ 0 H 1 − I
c c−
(3.4)
with I denoting the identity matrix. The entries of M are given by 1
P~ m μP~ n μμdμ; 1 m p δm;n−1 2 2 m − 12 m 1
The solution of boundary value problem (2.1) is discontinuous for μ 0 on the boundaries τ τ1 and τ τ2 . This discontinuity decays exponentially away from these boundaries due to scattering and absorption. In light of this discontinuity, we use the double spherical harmonics (DP N ) method to solve boundary value problem (2.1). The DP N method explicitly takes into account this discontinuity by defining I piecewise over the half-ranges, −1 ≤ μ < 0 and 0 < μ ≤ 1. We describe the DP N method in detail below. Let I μ; τ Iμ; τ and Q μ; τ Qμ; τ for 0 < μ ≤ 1, so that boundary value problem (2.1) becomes 1
(3.2)
where δmn denotes the Kronecker delta. We seek the solution of boundary value problem (3.1) as expansions of the form
M m1;n1
3. DOUBLE SPHERICAL HARMONICS METHOD
Z
P~ m μP~ n μdμ δmn ;
0
Z
where g is the anisotropy factor.
∂I ϖ μ I − 0 ∂τ 2
hμ; −μ0 I μ0 ; τ
I μ; τ1 f μ;
Consequently, I constant is a solution of Eq. (2.1a) when ϖ 0 1. In addition, 1 2
0
(3.1b)
We begin with the following boundary value problem for the one-dimensional RTE: ∂I ϖ μ I − 0 ∂τ 2
1
hμ; μ0 I − μ0 ; τdμ0 Q− μ; τ;
2. RADIATIVE TRANSPORT EQUATION
Z
Z
629
0
m1 δmn p δm;n1 ; 2 m 12 m 3
(3.5)
for m 0; 1; …, and n 0; 1; …, where we have made use of the recursion relation for Legendre polynomials: n 1P n1 x 2n 1xP n x − nP n−1 x, for −1 ≤ x ≤ 1. The matrix M is tridiagonal and symmetric. The entries of H 1 are given by H 1 m1;n1
1 2
Z
1 0
P~ m μ
Z
1
0
hμ; μ0 P~ n μ0 dμ0 dμ;
(3.6)
and the entries of H 2 are given by hμ; μ0 I μ0 ; τ
0 hμ; −μ0 I − μ0 ; τdμ0
Q μ; τ;
(3.1a)
H 2 m1;n1
1 2
Z
1 0
P~ m μ
Z 0
1
hμ; −μ0 P~ n μ0 dμ0 dμ;
(3.7)
630
J. Opt. Soc. Am. A / Vol. 31, No. 3 / March 2014
C. Sandoval and A. D. Kim
for m 0; 1; …, and n 0; 1; …. The matrices H 1 and H 2 are also symmetric. The components of the vectors q are q q 1 ; q2 ; …, where qn are defined as q n1 τ
Z
1 0
Q μ; τP~ n μdμ;
(3.8)
for n 0; 1; …. Boundary conditions (3.1c) and (3.1d) are prescribed over the half-range, 0 < μ ≤ 1. Consequently, for the DP N method, we prescribe the boundary conditions c τ1 f ;
(3.9a)
c− τ2 f − ;
(3.9b)
F τ
f −n1
Z
1
0
Z 0
1
f μP~ n μdμ;
(3.10)
F
M λ 0
0 −M
u ϖ 0 H 1 − I v ϖ 0 H 2
ϖ 0 H 2 ϖ 0 H 1 − I
u ; v
(3.12)
resulting from substituting c eλτ u and c− eλτ v into Eq. (3.4) with q 0. The eigenvalues, λ, and eigenvectors, u; v, are real. Here, we have introduced the notation u; v uT ; vT T , where the superscript T denotes the transpose. If λ is an eigenvalue with eigenvector u; v, then −λ is an eigenvalue with eigenvector v; u. In light of this symmetry, we index and order the eigenvalues according to λ−N1 ≤ ≤ λ−1 < 0 < λ1 ≤ ≤ λN1 ;
(3.13)
where λ−j −λj . We denote the eigenvectors corresponding to λj as uj ; vj for j 1; 2; . When ϖ 0 1, λ 0 is an eigenvalue, but otherwise λ ≠ 0, as indicated in Eq. (3.13). We make use of these properties in the analysis that follows.
4. GENERALIZED KUBELKA–MUNK SYSTEM The key step in deriving KM theory from the RTE lies in analyzing the DP N system with N 1, which we call the DP 1 system. To see why this relationship is important for this derivation, recall that KM theory provides the governing equations for the positive- and negative-going fluxes defined as
1 1 c0 p c 1 ; 2 3
1 1 p c0 c1 : G 2 3
(4.1)
(4.2)
(4.3)
Then, we have the following linear transformations
(3.11)
respectively, for n 0; 1; …. The system given by Eq. (3.4), subject to boundary conditions (3.9), is infinite dimensional. We compute an approximate solution by truncating the series in Eq. (3.3) after n N and neglecting c n for n > N. In doing so, the system given by Eq. (3.4) has dimension 2N 1. This finite dimensional system, which we call the DP N system, is suitable for numerical solution. Suppose we have performed this truncation. Consider the following 2N 1 × 2N 1 generalized eigenvalue problem:
I μ; τμdμ:
respectively. Thus, we require knowledge of only c 0 and c1 to compute F . The smallest DP N system involving just c0 and c 1 is the DP 1 system. Let us introduce the auxiliary quantities G defined as
f − −μP~ n μdμ;
1
0
Substituting Eq. (3.3) pinto Eq. (4.1), recognizing that μ 1∕2P~ 0 μ P~ 1 μ∕ 3, and making use of the orthogonality relation (3.2), we find that F in terms of the expansion coefficients c n are given by
where f and f − are the vectors whose components are given by f n1
Z
F G
"
1 2 1 p 2 3
1 p 2 3 1 2
#
c 0 : c 1
(4.4)
The 2 × 2 matrix in Eq. (4.4) is the matrix M whose entries are defined in Eq. (3.5) for m; n 0; 1. Let y F G T . Using the inverse transformation c M −1 y , we transform the DP 1 system given in Eq. (3.4) to
S1 I 0 d y S2 0 −I dτ y−
S2 S1
y q ; y− q−
(4.5)
with S 1 ϖ 0 H 1 M −1 − M −1 and S 2 ϖ 0 H 2 M −1 . By leftmultiplying boundary conditions (3.9) by M, we obtain the boundary conditions y τ1 Mf ;
(4.6a)
y− τ2 Mf − ;
(4.6b)
We call Eq. (4.5), the gKM equations, because it is a system governing the fluxes, F , and the auxiliary quantities, G . Boundary conditions for the gKM equations are given in Eq. (4.6). The gKM equations are valid when the truncation associated with the DP 1 system is valid. Since the DP N approximation is based on Eq. (3.3), which represents I as an expansion in orthogonal polynomials, the DP 1 system gives the optimal, piecewise linear representation of the specific intensity in the least-squares sense. Consider the 4 × 4 eigenvalue problem I 0 u~ S1 λ~ 0 −I v~ S2
S2 S1
u~ ; v~
(4.7)
associated to homogeneous gKM equations. Just as with the eigenvalue problem given in Eq. (3.12), if λ~ is an eigenvalue ~ v~ , then −λ~ is an eigenvalue with with eigenvector u; ~ Thus, we index and order the eigenvalues eigenvector ~v; u. according to
C. Sandoval and A. D. Kim
Vol. 31, No. 3 / March 2014 / J. Opt. Soc. Am. A
λ~ −2 ≤ λ~ −1 < 0 < λ~ 1 ≤ λ~ 2 ;
(4.8)
with λ~ −j −λ~ j . The eigenvectors corresponding to λ~ j are denoted by u~ j ; v~ j for j 1; 2. Consequently, the eigenvectors for λ~ −j are given by ~vj ; u~ j for j 1; 2. Thus, the general form for the solution of Eq. (4.5) is given by
y y−
2 X u~ j v~ a~ j τ j b~ j τ ; v~ j u~ j
(4.9)
termined by imposing that the solution satisfies the boundary conditions.
5. DERIVING THE KUBELKA–MUNK EQUATIONS The gKM equations are a 4 × 4 system of first-order differential equations. Its general solution given in Eq. (4.9) is a superposition of four linearly independent solutions. From the homogeneous solution, we identify two scales: λ~ 1 and λ~ 2 . In contrast, the KM equations given in Eq. (1.1) is a first-order p system with only one scale: KK 2S. To derive the KM equations, we study the DP 1 system for the case of strong multiple scattering. When multiple scattering is dominant, ϖ 0 ∼ 1. For that case, we compute an asymptotic solution of the generalized eigenvalue problem (3.12) for the DP 1 system from which we derive the KM equations. To make the assumption of strong multiple scattering precise, we introduce the small, dimensionless parameter, 0 < ϵ ≪ 1, and write ϖ 0 1 − ϵ2 , where ϵ2 μa ∕μs μa . For that case, after rearranging terms, Eq. (3.12) becomes
1 u H 2 u H 2 ϵ v v H 2 H 1 1 2 u H −I H : H 2 H 1 − I v
M λ 0
0 −M
(5.1)
We obtain below an asymptotic solution to this generalized eigenvalue problem in the limit as ϵ → 0 . When ϵ 0, Eq. (5.1) reduces to
M λ 0
now use perturbation theory to obtain higher order corrections to the λ 0 eigenvalue. Because of the symmetry property given in Eq. (3.13), we will therefore find the smallest (in magnitude) eigenvalues, λ1 , and their corresponding eigenvectors, given by u1 ; v1 and v1 ; u1 , respectively. We seek these higher order corrections in the form λ1 ϵλ01 ϵ2 λ001 Oϵ3 and u;v ˆe1 ; eˆ 1 ϵu0 1 ;v0 1 ϵ2 u00 1 ;v00 1 Oϵ3 . Substituting these asymptotic expansions into Eq. (5.1), and collecting like-powers of ϵ, we find to Oϵ that
j1
where the functions aj τ and bj τ for j 1; 2 are to be determined. For the homogeneous problem, in which q 0, ~ these functions are given by a~ j τ eλj τ−τ2 α~ j , and b~ j τ ~ e−λj τ−τ1 β~ j where the coefficients α~ j and β~ j for j 1; 2 are de-
1 u 0 H −I v −M H 2
u H 2 : H 1 − I v
(5.2)
We recognize that this is the case when ϖ 0 1 for Eq. (2.1a), for which I constant is a solution. This constant solution corresponds to the solution λ 0 and u v eˆ 1 , where eˆ 1 1; 0T . Substituting this solution into Eq. (5.2) yields H 1 eˆ 1 H 2 eˆ 1 eˆ 1 . In other words, the sum of the first columns of H 1 and H 2 are equal to eˆ 1 , which is a consequence of the normalization given in Eq. (2.2). The eigenvalue λ 0 is a double root of the characteristic polynomial associated with Eq. (5.2). By factoring out this double root from the characteristic polynomial, we obtain a quadratic equation with roots λ2 , where λ2 O1. We
631
λ01
M eˆ 1 −M eˆ 1
H 1 − I H 2
H 2 H 1 − I
u01 : 0 v1
(5.3)
For Eq. (5.3) to have a unique solution, its left-hand side must be orthogonal to ˆe1 ; eˆ 1 . Multiplying the left-hand side of Eq. (5.3) by ˆeT1 ; eˆ T1 , we obtain zero identically, so that this condition is satisfied. Next, we recognize that m1 M eˆ 1 is the first column of M, whose entries, according to Eq. (3.5), correspond to the projection of μ onto P~ m μ for m 0; 1. Therefore, the vector m1 ; −m1 corresponds to expansion coefficients for the functions I μ. According to Eq. (2.3), μ (over the whole range −1; 1) is an eigenfunction of the integral operator in Eq. (2.1a), with eigenvalue g. Consequently, the vector m1 ; −m1 is an eigenvector of the matrix in the right-hand side of Eq. (5.3), with eigenvalue −1 − g. Therefore, we find that
u01 v01
−
λ01 m1 : 1 − g −m1
(5.4)
To determine λ01 , we continue on to Oϵ2 , in which we find that " λ001
M
0
#"
eˆ 1
#
" λ01
M
0
#
u01
v01 eˆ 1 0 −M 0 −M # 00 # " " u1 eˆ 1 H 1 − I H 2 H 1 H 2 : (5.5) v001 eˆ 1 H 2 H 1 H 2 H 1 − I
The left-hand side of Eq. (5.5) must be orthogonal to ˆe1 ; eˆ 1 . By left-multiplying ˆeT1 ; eˆ T1 to the left-hand side of Eq. (5.5), setting that result to zero, substituting Eq. (5.4), and making use of the fact that H 1 eˆ 1 H 2 eˆ 1 eˆ 1 , we obtain, upon solving for λ02 1 , that λ02 1
1 − g : ∥m1 ∥2
(5.6)
p According to Eq. (3.5), m1 1∕2; 1∕2 3T . Thus, ∥m1 ∥2 p 1∕3, and so λ01 31 − g. p Thus, we have determined that λ1 ϵ 31 − g Oϵ2 is an eigenvalue with corresponding eigenvector u1 ; v1 p p ˆe1 − ϵ 3∕1 − gm1 ; eˆ 1 ϵ 3∕1 − gm1 Oϵ2 . By the symmetry property of the eigenvalues given in Section 3, −λ1 p −ϵ 31 − g is also an eigenvalue with corresponding eigenp p vector v1 ;u1 ˆe1 ϵ 3∕1−gm1 ; eˆ 1 −ϵ 3∕1−gm1 Oϵ2 . Since λ1 Oϵ is much smaller than λ2 O1, the leading order behavior of c governed by Eq. (3.4) when q 0 is given by
632
J. Opt. Soc. Am. A / Vol. 31, No. 3 / March 2014
c u v ∼ 1 eλ1 τ−τ2 α1 1 e−λ1 τ−τ1 β1 − v1 u1 c
C. Sandoval and A. D. Kim
(5.7)
in τ1 < τ < τ2 , where α1 and β1 are undetermined, scalar coefficients. Using Eq. (4.2) and our asymptotic results for u1 and v1 , we find that
F F−
u¯ v¯ v¯ u¯
eλ1 τ−τ2 0
0
e−λ1 τ−τ1
α1 ; β1
(5.8)
p p with u¯ 1∕2 − ϵ∕ 31 − g and v¯ 1∕2 ϵ∕ 31 − g. Equation (5.8) is the general solution of the first-order system d F u¯ − v¯ dτ F
v¯ u¯
λ1 0
0 −λ1
u¯ v¯ v¯ u¯
−1
F − : F
6. NUMERICAL EXAMPLES (5.9)
p p Substituting λ1 ϵ 31 − g, ϵ 1 − ϖ 0 , and carrying out the matrix products in Eq. (5.9), we obtain the KM equations, given in Eq. (1.1), with K 21 − ϖ 0 ;
(5.10)
3 S 1 − g − 1 − ϖ 0 : 4
(5.11)
Let K~ μs μa K and S~ μs μa S. Using ϖ 0 μs ∕μs μa and Eq. (5.10), we find that K~ 2μa . Similarly, using Eq. (5.11), we find that S~ 3∕4μs 1 − g − 1∕4μa 1–3g. The physically correct values for K~ and S~ have been the subject of several studies (see [8,11,14,16] and references contained therein). Thennadil [16] summarizes different values found in the literature by introducing the parameters x and y and a general form for these coefficients: K~ 2μa and S~ yμs 1 − g − xμa . For example, Gate [7] gives x 0, and y 3∕4. Brinkworth [6] gives x 1 and y 3∕4. Star et al. [9] give x 1∕4 and y 3∕4, which resembles the results above closest. The difference in the value for S~ obtained above in that it has an additional term that is proportional to the product μa g. For spectroscopy applications, KM theory is useful because one can compute an explicit expression for the diffuse reflectance from a slab of infinite thickness denoted by R∞ . The KM transform, 1 − R∞ 2 K~ ; ~S 2R∞
that the leading order asymptotic behavior given in Eq. (5.7) is defined on the open interval τ1 < τ < τ2 . This leading order asymptotic behavior is equivalent to neglecting terms with λ~ 2 in the solution of the gKM equations. Consequently, we have neglected two of the four linearly independent solutions for F in gKM to derive KM. It is for this reason that one must introduce the so-called four-flux theory to deal with collimated sources [2], for example. The addition of those two fluxes in four-flux theory is compensating for neglecting those other two solutions that are in the gKM equations.
(5.12)
~ and S~ [16]. We have gives a rational expression relating R∞ , K, ~ ~ determined K and S in terms of the fundamental optical properties: μa , μs , and g above. Thus, we have established a connection between R∞ and the fundamental optical properties. Moreover, because one can compute the reflectance and transmittance for a finite slab from R∞ [8], both of which are measurable quantities, one can readily extend these results to determine the fundamental optical properties of the slab. Thus, our results provide a systematic improvement to the application of KM theory to determine optical properties of a finite slab from measurable quantities. This derivation of the KM equations is based on the assumption of strong multiple scattering. For that case, we find a separation of scales corresponding to λ1 ≪ λ2 . Notice
In this section, we compare results from numerical solutions of the RTE, the gKM equations, and the KM equations. We first study numerical solutions of these equations for a homogeneous problem. Then, we study the nonhomogeneous problem for the diffuse flux. A. Homogeneous Problem We begin by studying the following homogeneous boundary value problem for the RTE μ
∂I ϖ I 0 ∂τ 2
Z
1
−1
hμ; μ0 Iμ0 ; τdμ0 2
Iμ; τ1 e−100μ−1
in τ1 < τ < τ2 ;
on 0 < μ ≤ 1;
Iμ; τ2 1 on − 1 ≤ μ < 0:
(6.1a)
(6.1b) (6.1c)
To solve boundary value problem (6.1), we use the discrete ordinate method described in Chapter 11 of [2]. In particular, we use the Henyey–Greenstein scattering phase function to compute the redistribution function defined as 21 − g2 2b p E ; (6.2) ab πa − b a b pp where a 1 g2 − 2gμμ0 , b 2g 1 − μ2 1 − μ02 , and Ek is the complete elliptic integral of the second kind. Upon numerical solution of boundary value problem (6.1), we compute F defined in Eq. (4.1) using the same Gauss-Legendre quadrature rule used for the discrete ordinate method. For the gKM system, we calculated numerically the entries of H 1 and H 2 using h defined in Eq. (6.2) and GaussLegendre quadrature. We also compute the vectors f defined in Eqs. (3.10) and (3.11) using Gauss-Legendre quadrature. With those matrices and vectors defined, we solve numerically the generalized eigenvalue problem (4.7) and obtain λ~ j , u~ j , and v~ j for j 1; 2. The solution of the gKM equations ~ is given by Eq. (4.9) with a~ j τ eλj τ−τ2 α~ j , and b~ j τ hμ; μ0
e−λj τ−τ1 β~ j . Evaluating this solution at τ τ1 and substituting that result into Eq. (4.6a), we find that ~
2 n X
o ~ u~ j e−λj τ2 −τ1 α~ j v~ j β~ j Mf :
(6.3)
j1
Evaluating the solution at τ τ2 and substituting that result into Eq. (4.6b), we find that
C. Sandoval and A. D. Kim 2 n X
Vol. 31, No. 3 / March 2014 / J. Opt. Soc. Am. A
o ~ v~ j α~ j u~ j e−λj τ2 −τ1 β~ j Mf − :
(6.4)
j1
Combining Eqs. (6.3) and (6.4) yields a 4 × 4 linear system of equations for α~ j and β~ j for j 1; 2. With α~ j and β~ j determined, we compute F through evaluation of F τ
2 n o X ~ ~ U~ 1j eλj τ−τ2 α~ j V~ 1j e−λj τ−τ1 β~ j ;
(6.5)
j1
F − τ
2 n X
o ~ ~ V~ 1j eλj τ−τ2 α~ j U~ 1j e−λj τ−τ1 β~ j :
(6.6)
j1
Here, U~ 1j denotes the first entry of u~ j , and V~ 1j denotes the first entry of v~ j . We solve the KM system (1.1) using K and S defined in Eqs. (5.10) and (5.11), respectively. We supplement this system with the boundary conditions Z1 F τ1 f μμdμ; (6.7) 0
F − τ2
Z
1
0
f − −μμdμ:
(6.8)
633
We compute the integrals in Eqs. (6.7) and (6.8) using GaussLegendre quadrature. In Fig. 1 we show results comparing solutions to (6.1) computed using the RTE (solid curves), the gKM (dashed curves), and the KM equations (dot-dashed curves). Here, we have set τ1 0, τ2 1, ϖ 0 0.99, and g 0. The top plot shows F as a function of τ, and the bottom plot shows F − as a function of τ. The gKM theory gives a more accurate solution compared to KM theory. It is nearly indistinguishable from that computed through the solution of the RTE. In Fig. 2 we plot the results when ϖ 0 0.99 and g 0.8. For this case, scattering is forward-peaked. Consequently, the specific intensity is not approximated well by a piecewise linear approximation. Thus, the DP 1 is not as good of an approximation, which, in turn, means that the gM and KM equations are less accurate. This is evidenced by the fact that both the gKM and KM results are less accurate than those shown in Fig. 1. Nonetheless, Fig. 2 shows that the gKM equations are still much more accurate than the KM equations. Next, we investigate reducing the value of ϖ 0 . For this case, the KM equations are known to be inaccurate [8,16,18]. Figure 3 shows results for ϖ 0 0.5 and g 0. Although scattering is isotropic, so that g 0, we find that the KM equations are much less accurate than the gKM equations. The error is worse for F τ than for F − τ presumably because it must also take into account the nonconstant boundary condition
0.3 RTE gKM KM
0.25
0.18 RTE gKM KM
0.16
F+(τ)
0.2
+
F (τ)
0.14
0.15
0.12
0.1 0.1
0.05
0
0.2
0.4
0.6
0.8
1
0.08
τ
0
0.2
0.4
0.55 RTE gKM KM
0.5
τ
0.6
0.8
1
0.6
0.8
1
0.52 RTE gKM KM
0.5
0.45
−
F (τ)
−
F (τ)
0.48
0.4
0.46 0.44
0.35 0.42 0
0.2
0.4
0.6
0.8
1
τ
F τ
F − τ
Fig. 1. Comparisons of (top) and (bottom) computed from solutions of boundary value problem (6.1) using the radiative transport equation (RTE), the generalized Kubelka–Munk (gKM) equations, and the Kubelka–Munk (KM) equations. Here, τ1 0, τ2 1, ϖ 0 0.99, and g 0.
0.4
0
0.2
0.4
τ Fig. 2. Comparisons of F τ (top) and F − τ (bottom) computed from solutions of boundary value problem (6.1) using the RTE, the gKM equations, and the KM equations. Here, τ1 0, τ2 1, ϖ 0 0.99, and g 0.8.
634
J. Opt. Soc. Am. A / Vol. 31, No. 3 / March 2014
C. Sandoval and A. D. Kim
(6.1b). When scattering is also forward-peaked, as shown in Fig. 4, we find that the KM equations produce negative values for F , which are unphysical. However, even with ϖ 0 0.5, the gKM equations yield accurate results. Even at a lower albedo (ϖ ∼ 0), the gKM equations yield qualitatively correct results. Although we do not show these results here, we have found that computing solutions with smaller (τ2 − τ1 ∼ 0.1) and larger (τ2 − τ1 ∼ 10) optical thicknesses yield similar results to those shown here. B. Nonhomogeneous Problem for the Diffuse Flux We now study the boundary value problem μ
∂I ϖ I 0 ∂τ 2
Z
1
−1
hμ; μ0 Iμ0 ; τdμ0
in τ1 < τ < τ2 ;
(6.9a)
Iμ; τ1 δμ − 1 on 0 < μ ≤ 1;
(6.9b)
Iμ; τ2 0 on − 1 ≤ μ < 0:
(6.9c)
This boundary value problem governs a plane wave with unit incident flux normally on a plane-parallel slab. It is a fundamental problem in the study of radiative transfer [1,2]. A complication in solving boundary value problem (6.9) is the Dirac delta function appearing in boundary condition (6.9b). To study this problem, one represents the specific intensity as
the sum I I ri I d . Here, I ri denotes the reduced intensity satisfying μ
(6.10a)
I ri μ; τ1 δμ − 1 on 0 < μ ≤ 1;
(6.10b)
I ri μ; τ2 0 on − 1 ≤ μ < 0:
(6.10c)
The solution of boundary value problem (6.10) is readily found to be I ri δμ − 1e−τ−τ1 ∕μ . The diffuse intensity, I d , satisfies the following nonhomogeneous, boundary value problem for the RTE: μ
∂I d ϖ Id 0 ∂τ 2
Z
1
−1
hμ; μ0 I d μ0 ; τdμ0 Qri ;
(6.11a)
I d μ; τ1 0 on 0 < μ ≤ 1;
(6.11b)
I d μ; τ2 0 on − 1 ≤ μ < 0:
(6.11c)
The nonhomogeneous term, Qri 1∕2ϖ 0 hμ; 1e−τ−τ1 , in Eq. (6.11a) comes from inserting I ri into the integral operation in the RTE. 0.1
0.11 RTE gKM KM
0.1
0.08 0.06 0.04
+
+
F (τ)
F (τ)
0.09
0.02
0.08
0 −0.02
RTE gKM KM
0.07 −0.04 0.06
∂I ri I ri 0; ∂τ
−0.06 0
0.2
0.4
τ
0.6
0.8
0
1
0.2
0.4
τ
0.6
0.8
1
0.6
0.8
1
0.55
0.55 RTE gKM KM
0.5
RTE gKM KM
0.5
0.45
0.45 F (τ)
0.35
0.4
−
−
F (τ)
0.4
0.35
0.3
0.3
0.25
0.25
0.2
0.2 0
0.2
0.4
0.6
0.8
1
τ Fig. 3. Comparisons of F τ (top) and F − τ (bottom) computed from solutions of boundary value problem (6.1) using the RTE, the gKM equations, and the KM equations. Here, τ1 0, τ2 1, ϖ 0 0.5, and g 0.
0
0.2
0.4 τ
Fig. 4. Comparisons of F τ (top) and F − τ (bottom) computed from solutions of boundary value problem (6.1) using the RTE, the gKM equations, and the KM equations. Here, τ1 0, τ2 1, ϖ 0 0.5, and g 0.8.
C. Sandoval and A. D. Kim
Vol. 31, No. 3 / March 2014 / J. Opt. Soc. Am. A
results for y y H yP at τ τ1 and substituting that result into Eq. (4.6a) with f 0, we find that
We compute the diffuse fluxes, defined as F d τ
Z
1 0
I d μ; τμdμ;
(6.12)
by solving boundary value problem (6.11) directly and by solving the gKM for this problem. To solve boundary value problem (6.11), we use the same method that we use to solve boundary value problem (6.1). The gKM equations for this problem take the form
I 0 d y S S 1 2 0 −I dτ y− S2 S1
ϖ 0 h −τ−τ1 y e ; (6.13) y− 2 h−
where the entries of the vectors h are given by Z1 h hμ; 1P n μdμ; n 0; 1: n1
−I − S 1 −S 2
−S 2 I − S1
η η−
ϖ h 0 − : 2 h
2 n X
o ~ v~ j α~ j u~ j e−λj τ2 −τ1 β~ j −η− e−τ2 −τ1 :
(6.17)
j1
Combining Eqs. (6.3) and (6.4) yields a 4 × 4 linear system of equations for α~ j and β~ j for j 1; 2. With α~ j and β~ j determined, we compute F through evaluation of 2 n X
o ~ ~ −τ−τ1 ; (6.18) U~ 1j eλj τ−τ2 α~ j V~ 1j e−λj τ−τ1 β~ j η 1 e
j1
F − τ
2 n o X ~ ~ V~ 1j eλj τ−τ2 α~ j U~ 1j e−λj τ−τ1 β~ j η−1 e−τ−τ1 ;
(6.19) with η 1 denoting the first entries of η . In Fig. 5, we compare the diffuse fluxes computed from solutions of (6.11) by solving the RTE (solid curves) and
0.05
0.1
0.04
0.08
0.03
0.06
F (τ)
0.02
+
+
Evaluating the solution at τ τ2 and substituting that result into Eq. (4.6b) with f − 0, we find that
(6.15)
The homogeneous solution is given by Eq. (4.9) with ~ ~ a~ j τ eλj τ−τ2 α~ j , and b~ j τ e−λj τ−τ1 β~ j . Evaluating our
0.01
0.04 0.02
0
0
RTE gKM 0
0.2
0.4
τ
0.6
0.8
−0.02
1
0.06
0.016
0.05
0.014
RTE gKM 0
0.2
0.4
τ
0.6
0.8
1
0.6
0.8
1
0.012
0.04
0.01
F (τ)
0.03
−
−
F (τ)
(6.16)
j1
y H
−0.01
o ~ u~ j e−λj τ2 −τ1 α~ j v~ j β~ j −η :
j1
F τ
We compute the solution of Eq. (6.13) as the sum of the homo geneous and particular solutions: y y H yP . We seek the −τ−τ 1 , where η particular solution in the form yP η e satisfies the 4 × 4 linear system
F (τ)
2 n X
(6.14)
0
635
0.02
0.008 0.006
0.01
0.004
0 −0.01
RTE gKM 0
0.2
0.002 0.4
0.6
0.8
1
τ Fig. 5. Comparisons of F τ (top) and F − τ (bottom) computed from solutions of boundary value problem (6.11) using the RTE, and the gKM equations. Here, τ1 0, τ2 1, ϖ 0 0.99, and g 0.
0
RTE gKM 0
0.2
0.4
τ
Fig. 6. Comparisons of F τ (top) and F − τ (bottom) computed from solutions of boundary value problem (6.11) using the RTE and the gKM equations. Here, τ1 0, τ2 1, ϖ 0 0.99, and g 0.7.
636
J. Opt. Soc. Am. A / Vol. 31, No. 3 / March 2014
the gKM equations (dashed curves). Here, we have set τ1 0, τ2 1, ϖ 0 0.99, and g 0. Here, the gKM equations provide a very accurate approximation to the RTE. Figure 6 compares the diffuse fluxes for the case in which g 0.7. As expected, the approximation produced by the gKM equations are less accurate when scattering is forward-peaked than when it is isotropic. Nonetheless, the gKM equations still provide an accurate approximation to the RTE for the diffuse fluxes. Although we do not show these results here, we find that gKM equations still approximate the RTE in the manner shown in Figs. 5 and 6 with lower albedo values, and for both smaller and larger optical thicknesses.
7. CONCLUSION We have systematically derived the KM equations from the RTE. This derivation involves the application of the DP N method to solve the RTE. By transforming the DP 1 system of equations, we obtained the gKM equations. Then, by studying the DP 1 system of equations in the limit of strong multiple scattering, we derived the KM equations. For the KM equations, we have determined all coefficients explicitly. In particular, we have found that K~ 2μa and S~ 3∕4μs 1 − g − 1∕4μa 1–3g. These results agree closely with those found in the literature. However, the significant difference in the results found here is the term in S~ proportional to μa g. Moreover, we established this theory’s range of validity as a consequence of the asymptotic analysis used to derive it. The gKM equations can easily account for general boundary conditions and nonhomogeneous problems. The numerical results shown here demonstrate that the gKM equations are much better than the KM equations at approximating the solutions of the RTE. The complexity of the gKM equations is only slightly larger, since it is a 4 × 4 system rather than a 2 × 2 one. Since it is more accurate than the KM system of equations, more broadly applicable, and only slightly more difficult to solve, the gKM equations should be very useful.
REFERENCES 1. 2.
S. Chandrasekhar, Radiative Transfer (Dover, 1960). A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE, 1996).
C. Sandoval and A. D. Kim 3. P. Kubelka and F. Munk, “Ein beitrag zur optik der farbanstriche,” Z. Tech. Phys. 12, 593–601 (1931). 4. P. Kubelka, “New contributions to the optics of intensely lightscattering materials. Part I,” J. Opt. Soc. Am. 38, 448–457 (1948). 5. B. Philips-Invernizzi, D. Dupont, and C. Cazé, “Bibliographical review for reflectance of diffusing media,” Opt. Eng. 40, 1082–1092 (2001). 6. B. J. Brinkworth, “Interpretations of the Kubelka-Munk coefficients in reflection theory,” Appl. Opt. 11, 1434–1435 (1972). 7. L. F. Gate, “Comparison of the photon diffusion model and Kubelka-Munk equation with exact solution of the radiative transport equation,” Appl. Opt. 13, 236–238 (1974). 8. J. Nobbs, “Kubelka-Munk theory and the prediction of reflectance,” Rev. Prog. Coloration 15, 66–75 (1985). 9. W. M. Star, J. P. A. Marijnissen, and M. J. C. Van Gemert, “Light dosimetry in optical phantoms and in tissues: I. Multiple flux and transport theory,” Phys. Med. Biol. 33, 437–454 (1988). 10. W. E. Vargas and G. A. Niklasson, “Applicability conditions of the Kubelka-Munk theory,” Appl. Opt. 36, 5580–5586 (1997). 11. R. Molenaar, J. ten Bosch, and J. Zijp, “Determination of Kubelka-Munk scattering and absorption coefficients by diffuse illumination,” Appl. Opt. 38, 2068–2077 (1999). 12. L. Yang and B. Kruse, “Revised Kubelka-Munk theory. I. Theory and application,” J. Opt. Soc. Am. A 21, 1933–1941 (2004). 13. L. Yang, B. Kruse, and S. J. Miklavcic, “Revised Kubelka-Munk theory. II. Unified framework for homogeneous and inhomogeneous optical media,” J. Opt. Soc. Am. A 21, 1942–1952 (2004). 14. L. Yang and S. J. Miklavcic, “Revised Kubelka-Munk theory. III. A general theory of light propagation in scattering and absorptive media,” J. Opt. Soc. Am. A 22, 1866–1873 (2005). 15. P. Edström, “Examination of the revised Kubelka-Munk theory: considerations of modeling strategies,” J. Opt. Soc. Am. A 24, 548–556 (2007). 16. S. N. Thennadil, “Relationships between the Kubelka-Munk scattering and radiative transfer coefficients,” J. Opt. Soc. Am. A 25, 1480–1485 (2008). 17. M. Neuman and P. Edström, “Anisotropic reflectance from turbid media. I. Theory,” J. Opt. Soc. Am. A 27, 1032–1039 (2010). 18. M. L. Myrick, M. N. Simcock, M. Baranowski, H. Brooke, S. L. Morgan, and J. N. McCutcheon, “The Kubelka-Munk diffuse reflectance formula revisited,” Appl. Spectrosc. Rev. 46, 140–165 (2011). 19. B. Davison, Neutron Transport Theory (Oxford University, 1958). 20. K. M. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, 1967). 21. E. E. Lewis and W. F. Miller, Computational Methods of Neutron Transport (American Nuclear Society, 1993). 22. R. Aronson, “PN vs. double-PN approximations for highly anisotropic scattering,” Transp. Theory Stat. Phys. 15, 829–840 (1986).
J. Opt. Soc. Am. A / Vol. 31, No. 3 / March 2014
C. Sandoval and A. D. Kim
Deriving Kubelka–Munk theory from radiative transport Christopher Sandoval and Arnold D. Kim* Applied Mathematics Unit, School of Natural Sciences, University of California, Merced, 5200 North Lake Road, Merced, California 95343, USA *Corresponding author: [email protected] Received November 12, 2013; accepted January 6, 2014; posted January 16, 2014 (Doc. ID 201164); published February 21, 2014 We derive Kubelka–Munk (KM) theory systematically from the radiative transport equation (RTE) by analyzing the system of equations resulting from applying the double spherical harmonics method of order one and transforming that system into one governing the positive- and negative-going fluxes. Through this derivation, we establish the theoretical basis of KM theory, identify all parameters, and determine its range of validity. Moreover, we are able to generalize KM theory to take into account general boundary sources and nonhomogeneous terms, for example. The generalized Kubelka–Munk (gKM) equations are also much more accurate at approximating the solution of the RTE. We validate this theory through comparison with numerical solutions of the RTE. © 2014 Optical Society of America OCIS codes: (290.4210) Multiple scattering; (290.7050) Turbid media; (030.5620) Radiative transfer; (000.3860) Mathematical methods in physics. http://dx.doi.org/10.1364/JOSAA.31.000628
1. INTRODUCTION Multiple scattering of light in a turbid medium is governed by the theory of radiative transport [1,2]. This theory provides a complete description of absorption, scattering, and radiation of light. However, applying radiative transport theory is challenging for practical problems because exact solutions of the radiative transport equation (RTE) are known only for simple problems. Even numerical solutions are challenging due to the large number of variables involved. Kubelka and Munk [3,4] developed a phenomenological theory that describes the flow of power in positive- and negative-going directions, which we denote by F , respectively, through a plane parallel slab of a scattering medium. The Kubelka–Munk (KM) equations are dF K SF SF − ; dτ
(1.1a)
dF − K SF − SF ; dτ
(1.1b)
−
where τ is the optical depth, and K and S are related to absorption and scattering, respectively. This so-called two-flux theory is intuitive and much simpler than the RTE. Hence, it has been applied to several practical problems (see [5] and references contained therein). However, the theoretical basis of KM theory is not entirely established [2]. Several studies [6–18] have addressed the theoretical basis for KM theory and its applicability to interpret measured data. However, a systematic derivation of the KM equations from the RTE is still needed to determine all coefficients, to determine its range of validity, and to provide a framework to derive extensions that will allow for broader application to practical problems. We present here a systematic derivation of the KM equations from the RTE. This derivation comes from analyzing 1084-7529/14/030628-09$15.00/0
the system of equations resulting from applying the double spherical harmonics method of order one to a general boundary value problem for the one dimensional RTE. This theoretical framework also provides a means to generalize KM theory to a larger class of problems including general boundary conditions and nonhomogeneous terms. For example, we identify the fundamental reason why a four-flux, rather than a two-flux, theory is needed to take into account collimated sources. We validate all aspects of this theory through comparisons with numerical solutions of the RTE. A key step in this derivation is the application of the double spherical harmonics method [19–22] to solve the general boundary value problem for the one dimensional RTE. The double spherical harmonics method uses two different expansions in the angular variable to represent the specific intensity over two distinct hemispheres. Doing so also provides a natural representation for the solution which, in general, is discontinuous at directions tangent to the parallel boundary planes. The double spherical harmonics method provides a natural way to prescribe boundary conditions since they are prescribed over half ranges. It is this piecewise representation of the specific intensity that allows for the systematic derivation of the KM equations governing the positive- and negativegoing fluxes. The remainder of this paper is as follows. In Section 2, we present the general boundary value problem for the one dimensional RTE. In Section 3, we apply the double spherical harmonics method to solve this boundary value problem for the RTE. In Section 4, we transform the double spherical harmonics system of order one to a 4 × 4 system, governing the positive- and negative-going fluxes and two auxiliary functions, which we call the generalized Kubelka–Munk (gKM) equations. In Section 5, we reduce this generalized system to the 2 × 2 system of KM equations, given in Eq. (1.1), assuming that multiple scattering is strong. We give numerical © 2014 Optical Society of America
C. Sandoval and A. D. Kim
Vol. 31, No. 3 / March 2014 / J. Opt. Soc. Am. A
results of comparisons between results of this theory and those of direct numerical solution of the RTE in Section 6. In Section 7, we give our conclusions.
−μ
∂I − ϖ I− − 0 ∂τ 2
1
−1
hμ; μ0 Iμ0 ; τdμ0 Qμ; τ in τ1 < τ < τ2 ; (2.1a)
Iμ; τ1 f μ on 0 < μ ≤ 1;
(2.1b)
Iμ; τ2 f − μ on − 1 ≤ μ < 0:
(2.1c)
Here, Iμ; τ is the specific intensity which depends on μ cos θ, where θ is the polar angle, and τ is the optical depth. We denote the albedo by ϖ 0 , which is defined in terms of the absorption and scattering coefficients, μa and μs , respectively, as ϖ 0 μs ∕μs μa . The function Qμ; τ is a source, and the functions f give sources on the boundaries, τ τ1 and τ τ2 , respectively, for all directions pointing into the plane-parallel slab. The redistribution function, h, gives the fraction of light scattered in direction μ due to light incident in direction μ0 . It is normalized according to Z
1 2
1
−1
hμ; μ0 dμ0 1:
(2.2)
Z
1
−1
hμ; μ0 μ0 dμ0 gμ;
(2.3)
(3.1c)
I − μ; τ2 f − −μ;
(3.1d)
Here, we have assumed that scattering is spherically symmetric, so that h−μ; −μ0 hμ; μ0 and hμ; −μ0 h−μ; μ0 . We p introduce the functions P~ n μ 2n 1P n 2μ − 1 defined over the half-range, 0 < μ ≤ 1, where fP n g is the set of Legendre polynomials. We call fP~ n g the set of normalized, half-range Legendre polynomials. They satisfy the following orthogonality relation: Z
1
I μ; τ
∞ X
~ c n τP n μ:
(3.3)
n0
Let c denote the vectors whose components are c T c 0 ; c1 ; . By substituting Eq. (3.3) and making use of orthogonality relation (3.2), we find that c satisfies the system d M dτ 0
0 −M
c c−
ϖ 0 H 1 − I ϖ 0 H 2 q ; q−
ϖ 0 H 2 ϖ 0 H 1 − I
c c−
(3.4)
with I denoting the identity matrix. The entries of M are given by 1
P~ m μP~ n μμdμ; 1 m p δm;n−1 2 2 m − 12 m 1
The solution of boundary value problem (2.1) is discontinuous for μ 0 on the boundaries τ τ1 and τ τ2 . This discontinuity decays exponentially away from these boundaries due to scattering and absorption. In light of this discontinuity, we use the double spherical harmonics (DP N ) method to solve boundary value problem (2.1). The DP N method explicitly takes into account this discontinuity by defining I piecewise over the half-ranges, −1 ≤ μ < 0 and 0 < μ ≤ 1. We describe the DP N method in detail below. Let I μ; τ Iμ; τ and Q μ; τ Qμ; τ for 0 < μ ≤ 1, so that boundary value problem (2.1) becomes 1
(3.2)
where δmn denotes the Kronecker delta. We seek the solution of boundary value problem (3.1) as expansions of the form
M m1;n1
3. DOUBLE SPHERICAL HARMONICS METHOD
Z
P~ m μP~ n μdμ δmn ;
0
Z
where g is the anisotropy factor.
∂I ϖ μ I − 0 ∂τ 2
hμ; −μ0 I μ0 ; τ
I μ; τ1 f μ;
Consequently, I constant is a solution of Eq. (2.1a) when ϖ 0 1. In addition, 1 2
0
(3.1b)
We begin with the following boundary value problem for the one-dimensional RTE: ∂I ϖ μ I − 0 ∂τ 2
1
hμ; μ0 I − μ0 ; τdμ0 Q− μ; τ;
2. RADIATIVE TRANSPORT EQUATION
Z
Z
629
0
m1 δmn p δm;n1 ; 2 m 12 m 3
(3.5)
for m 0; 1; …, and n 0; 1; …, where we have made use of the recursion relation for Legendre polynomials: n 1P n1 x 2n 1xP n x − nP n−1 x, for −1 ≤ x ≤ 1. The matrix M is tridiagonal and symmetric. The entries of H 1 are given by H 1 m1;n1
1 2
Z
1 0
P~ m μ
Z
1
0
hμ; μ0 P~ n μ0 dμ0 dμ;
(3.6)
and the entries of H 2 are given by hμ; μ0 I μ0 ; τ
0 hμ; −μ0 I − μ0 ; τdμ0
Q μ; τ;
(3.1a)
H 2 m1;n1
1 2
Z
1 0
P~ m μ
Z 0
1
hμ; −μ0 P~ n μ0 dμ0 dμ;
(3.7)
630
J. Opt. Soc. Am. A / Vol. 31, No. 3 / March 2014
C. Sandoval and A. D. Kim
for m 0; 1; …, and n 0; 1; …. The matrices H 1 and H 2 are also symmetric. The components of the vectors q are q q 1 ; q2 ; …, where qn are defined as q n1 τ
Z
1 0
Q μ; τP~ n μdμ;
(3.8)
for n 0; 1; …. Boundary conditions (3.1c) and (3.1d) are prescribed over the half-range, 0 < μ ≤ 1. Consequently, for the DP N method, we prescribe the boundary conditions c τ1 f ;
(3.9a)
c− τ2 f − ;
(3.9b)
F τ
f −n1
Z
1
0
Z 0
1
f μP~ n μdμ;
(3.10)
F
M λ 0
0 −M
u ϖ 0 H 1 − I v ϖ 0 H 2
ϖ 0 H 2 ϖ 0 H 1 − I
u ; v
(3.12)
resulting from substituting c eλτ u and c− eλτ v into Eq. (3.4) with q 0. The eigenvalues, λ, and eigenvectors, u; v, are real. Here, we have introduced the notation u; v uT ; vT T , where the superscript T denotes the transpose. If λ is an eigenvalue with eigenvector u; v, then −λ is an eigenvalue with eigenvector v; u. In light of this symmetry, we index and order the eigenvalues according to λ−N1 ≤ ≤ λ−1 < 0 < λ1 ≤ ≤ λN1 ;
(3.13)
where λ−j −λj . We denote the eigenvectors corresponding to λj as uj ; vj for j 1; 2; . When ϖ 0 1, λ 0 is an eigenvalue, but otherwise λ ≠ 0, as indicated in Eq. (3.13). We make use of these properties in the analysis that follows.
4. GENERALIZED KUBELKA–MUNK SYSTEM The key step in deriving KM theory from the RTE lies in analyzing the DP N system with N 1, which we call the DP 1 system. To see why this relationship is important for this derivation, recall that KM theory provides the governing equations for the positive- and negative-going fluxes defined as
1 1 c0 p c 1 ; 2 3
1 1 p c0 c1 : G 2 3
(4.1)
(4.2)
(4.3)
Then, we have the following linear transformations
(3.11)
respectively, for n 0; 1; …. The system given by Eq. (3.4), subject to boundary conditions (3.9), is infinite dimensional. We compute an approximate solution by truncating the series in Eq. (3.3) after n N and neglecting c n for n > N. In doing so, the system given by Eq. (3.4) has dimension 2N 1. This finite dimensional system, which we call the DP N system, is suitable for numerical solution. Suppose we have performed this truncation. Consider the following 2N 1 × 2N 1 generalized eigenvalue problem:
I μ; τμdμ:
respectively. Thus, we require knowledge of only c 0 and c1 to compute F . The smallest DP N system involving just c0 and c 1 is the DP 1 system. Let us introduce the auxiliary quantities G defined as
f − −μP~ n μdμ;
1
0
Substituting Eq. (3.3) pinto Eq. (4.1), recognizing that μ 1∕2P~ 0 μ P~ 1 μ∕ 3, and making use of the orthogonality relation (3.2), we find that F in terms of the expansion coefficients c n are given by
where f and f − are the vectors whose components are given by f n1
Z
F G
"
1 2 1 p 2 3
1 p 2 3 1 2
#
c 0 : c 1
(4.4)
The 2 × 2 matrix in Eq. (4.4) is the matrix M whose entries are defined in Eq. (3.5) for m; n 0; 1. Let y F G T . Using the inverse transformation c M −1 y , we transform the DP 1 system given in Eq. (3.4) to
S1 I 0 d y S2 0 −I dτ y−
S2 S1
y q ; y− q−
(4.5)
with S 1 ϖ 0 H 1 M −1 − M −1 and S 2 ϖ 0 H 2 M −1 . By leftmultiplying boundary conditions (3.9) by M, we obtain the boundary conditions y τ1 Mf ;
(4.6a)
y− τ2 Mf − ;
(4.6b)
We call Eq. (4.5), the gKM equations, because it is a system governing the fluxes, F , and the auxiliary quantities, G . Boundary conditions for the gKM equations are given in Eq. (4.6). The gKM equations are valid when the truncation associated with the DP 1 system is valid. Since the DP N approximation is based on Eq. (3.3), which represents I as an expansion in orthogonal polynomials, the DP 1 system gives the optimal, piecewise linear representation of the specific intensity in the least-squares sense. Consider the 4 × 4 eigenvalue problem I 0 u~ S1 λ~ 0 −I v~ S2
S2 S1
u~ ; v~
(4.7)
associated to homogeneous gKM equations. Just as with the eigenvalue problem given in Eq. (3.12), if λ~ is an eigenvalue ~ v~ , then −λ~ is an eigenvalue with with eigenvector u; ~ Thus, we index and order the eigenvalues eigenvector ~v; u. according to
C. Sandoval and A. D. Kim
Vol. 31, No. 3 / March 2014 / J. Opt. Soc. Am. A
λ~ −2 ≤ λ~ −1 < 0 < λ~ 1 ≤ λ~ 2 ;
(4.8)
with λ~ −j −λ~ j . The eigenvectors corresponding to λ~ j are denoted by u~ j ; v~ j for j 1; 2. Consequently, the eigenvectors for λ~ −j are given by ~vj ; u~ j for j 1; 2. Thus, the general form for the solution of Eq. (4.5) is given by
y y−
2 X u~ j v~ a~ j τ j b~ j τ ; v~ j u~ j
(4.9)
termined by imposing that the solution satisfies the boundary conditions.
5. DERIVING THE KUBELKA–MUNK EQUATIONS The gKM equations are a 4 × 4 system of first-order differential equations. Its general solution given in Eq. (4.9) is a superposition of four linearly independent solutions. From the homogeneous solution, we identify two scales: λ~ 1 and λ~ 2 . In contrast, the KM equations given in Eq. (1.1) is a first-order p system with only one scale: KK 2S. To derive the KM equations, we study the DP 1 system for the case of strong multiple scattering. When multiple scattering is dominant, ϖ 0 ∼ 1. For that case, we compute an asymptotic solution of the generalized eigenvalue problem (3.12) for the DP 1 system from which we derive the KM equations. To make the assumption of strong multiple scattering precise, we introduce the small, dimensionless parameter, 0 < ϵ ≪ 1, and write ϖ 0 1 − ϵ2 , where ϵ2 μa ∕μs μa . For that case, after rearranging terms, Eq. (3.12) becomes
1 u H 2 u H 2 ϵ v v H 2 H 1 1 2 u H −I H : H 2 H 1 − I v
M λ 0
0 −M
(5.1)
We obtain below an asymptotic solution to this generalized eigenvalue problem in the limit as ϵ → 0 . When ϵ 0, Eq. (5.1) reduces to
M λ 0
now use perturbation theory to obtain higher order corrections to the λ 0 eigenvalue. Because of the symmetry property given in Eq. (3.13), we will therefore find the smallest (in magnitude) eigenvalues, λ1 , and their corresponding eigenvectors, given by u1 ; v1 and v1 ; u1 , respectively. We seek these higher order corrections in the form λ1 ϵλ01 ϵ2 λ001 Oϵ3 and u;v ˆe1 ; eˆ 1 ϵu0 1 ;v0 1 ϵ2 u00 1 ;v00 1 Oϵ3 . Substituting these asymptotic expansions into Eq. (5.1), and collecting like-powers of ϵ, we find to Oϵ that
j1
where the functions aj τ and bj τ for j 1; 2 are to be determined. For the homogeneous problem, in which q 0, ~ these functions are given by a~ j τ eλj τ−τ2 α~ j , and b~ j τ ~ e−λj τ−τ1 β~ j where the coefficients α~ j and β~ j for j 1; 2 are de-
1 u 0 H −I v −M H 2
u H 2 : H 1 − I v
(5.2)
We recognize that this is the case when ϖ 0 1 for Eq. (2.1a), for which I constant is a solution. This constant solution corresponds to the solution λ 0 and u v eˆ 1 , where eˆ 1 1; 0T . Substituting this solution into Eq. (5.2) yields H 1 eˆ 1 H 2 eˆ 1 eˆ 1 . In other words, the sum of the first columns of H 1 and H 2 are equal to eˆ 1 , which is a consequence of the normalization given in Eq. (2.2). The eigenvalue λ 0 is a double root of the characteristic polynomial associated with Eq. (5.2). By factoring out this double root from the characteristic polynomial, we obtain a quadratic equation with roots λ2 , where λ2 O1. We
631
λ01
M eˆ 1 −M eˆ 1
H 1 − I H 2
H 2 H 1 − I
u01 : 0 v1
(5.3)
For Eq. (5.3) to have a unique solution, its left-hand side must be orthogonal to ˆe1 ; eˆ 1 . Multiplying the left-hand side of Eq. (5.3) by ˆeT1 ; eˆ T1 , we obtain zero identically, so that this condition is satisfied. Next, we recognize that m1 M eˆ 1 is the first column of M, whose entries, according to Eq. (3.5), correspond to the projection of μ onto P~ m μ for m 0; 1. Therefore, the vector m1 ; −m1 corresponds to expansion coefficients for the functions I μ. According to Eq. (2.3), μ (over the whole range −1; 1) is an eigenfunction of the integral operator in Eq. (2.1a), with eigenvalue g. Consequently, the vector m1 ; −m1 is an eigenvector of the matrix in the right-hand side of Eq. (5.3), with eigenvalue −1 − g. Therefore, we find that
u01 v01
−
λ01 m1 : 1 − g −m1
(5.4)
To determine λ01 , we continue on to Oϵ2 , in which we find that " λ001
M
0
#"
eˆ 1
#
" λ01
M
0
#
u01
v01 eˆ 1 0 −M 0 −M # 00 # " " u1 eˆ 1 H 1 − I H 2 H 1 H 2 : (5.5) v001 eˆ 1 H 2 H 1 H 2 H 1 − I
The left-hand side of Eq. (5.5) must be orthogonal to ˆe1 ; eˆ 1 . By left-multiplying ˆeT1 ; eˆ T1 to the left-hand side of Eq. (5.5), setting that result to zero, substituting Eq. (5.4), and making use of the fact that H 1 eˆ 1 H 2 eˆ 1 eˆ 1 , we obtain, upon solving for λ02 1 , that λ02 1
1 − g : ∥m1 ∥2
(5.6)
p According to Eq. (3.5), m1 1∕2; 1∕2 3T . Thus, ∥m1 ∥2 p 1∕3, and so λ01 31 − g. p Thus, we have determined that λ1 ϵ 31 − g Oϵ2 is an eigenvalue with corresponding eigenvector u1 ; v1 p p ˆe1 − ϵ 3∕1 − gm1 ; eˆ 1 ϵ 3∕1 − gm1 Oϵ2 . By the symmetry property of the eigenvalues given in Section 3, −λ1 p −ϵ 31 − g is also an eigenvalue with corresponding eigenp p vector v1 ;u1 ˆe1 ϵ 3∕1−gm1 ; eˆ 1 −ϵ 3∕1−gm1 Oϵ2 . Since λ1 Oϵ is much smaller than λ2 O1, the leading order behavior of c governed by Eq. (3.4) when q 0 is given by
632
J. Opt. Soc. Am. A / Vol. 31, No. 3 / March 2014
c u v ∼ 1 eλ1 τ−τ2 α1 1 e−λ1 τ−τ1 β1 − v1 u1 c
C. Sandoval and A. D. Kim
(5.7)
in τ1 < τ < τ2 , where α1 and β1 are undetermined, scalar coefficients. Using Eq. (4.2) and our asymptotic results for u1 and v1 , we find that
F F−
u¯ v¯ v¯ u¯
eλ1 τ−τ2 0
0
e−λ1 τ−τ1
α1 ; β1
(5.8)
p p with u¯ 1∕2 − ϵ∕ 31 − g and v¯ 1∕2 ϵ∕ 31 − g. Equation (5.8) is the general solution of the first-order system d F u¯ − v¯ dτ F
v¯ u¯
λ1 0
0 −λ1
u¯ v¯ v¯ u¯
−1
F − : F
6. NUMERICAL EXAMPLES (5.9)
p p Substituting λ1 ϵ 31 − g, ϵ 1 − ϖ 0 , and carrying out the matrix products in Eq. (5.9), we obtain the KM equations, given in Eq. (1.1), with K 21 − ϖ 0 ;
(5.10)
3 S 1 − g − 1 − ϖ 0 : 4
(5.11)
Let K~ μs μa K and S~ μs μa S. Using ϖ 0 μs ∕μs μa and Eq. (5.10), we find that K~ 2μa . Similarly, using Eq. (5.11), we find that S~ 3∕4μs 1 − g − 1∕4μa 1–3g. The physically correct values for K~ and S~ have been the subject of several studies (see [8,11,14,16] and references contained therein). Thennadil [16] summarizes different values found in the literature by introducing the parameters x and y and a general form for these coefficients: K~ 2μa and S~ yμs 1 − g − xμa . For example, Gate [7] gives x 0, and y 3∕4. Brinkworth [6] gives x 1 and y 3∕4. Star et al. [9] give x 1∕4 and y 3∕4, which resembles the results above closest. The difference in the value for S~ obtained above in that it has an additional term that is proportional to the product μa g. For spectroscopy applications, KM theory is useful because one can compute an explicit expression for the diffuse reflectance from a slab of infinite thickness denoted by R∞ . The KM transform, 1 − R∞ 2 K~ ; ~S 2R∞
that the leading order asymptotic behavior given in Eq. (5.7) is defined on the open interval τ1 < τ < τ2 . This leading order asymptotic behavior is equivalent to neglecting terms with λ~ 2 in the solution of the gKM equations. Consequently, we have neglected two of the four linearly independent solutions for F in gKM to derive KM. It is for this reason that one must introduce the so-called four-flux theory to deal with collimated sources [2], for example. The addition of those two fluxes in four-flux theory is compensating for neglecting those other two solutions that are in the gKM equations.
(5.12)
~ and S~ [16]. We have gives a rational expression relating R∞ , K, ~ ~ determined K and S in terms of the fundamental optical properties: μa , μs , and g above. Thus, we have established a connection between R∞ and the fundamental optical properties. Moreover, because one can compute the reflectance and transmittance for a finite slab from R∞ [8], both of which are measurable quantities, one can readily extend these results to determine the fundamental optical properties of the slab. Thus, our results provide a systematic improvement to the application of KM theory to determine optical properties of a finite slab from measurable quantities. This derivation of the KM equations is based on the assumption of strong multiple scattering. For that case, we find a separation of scales corresponding to λ1 ≪ λ2 . Notice
In this section, we compare results from numerical solutions of the RTE, the gKM equations, and the KM equations. We first study numerical solutions of these equations for a homogeneous problem. Then, we study the nonhomogeneous problem for the diffuse flux. A. Homogeneous Problem We begin by studying the following homogeneous boundary value problem for the RTE μ
∂I ϖ I 0 ∂τ 2
Z
1
−1
hμ; μ0 Iμ0 ; τdμ0 2
Iμ; τ1 e−100μ−1
in τ1 < τ < τ2 ;
on 0 < μ ≤ 1;
Iμ; τ2 1 on − 1 ≤ μ < 0:
(6.1a)
(6.1b) (6.1c)
To solve boundary value problem (6.1), we use the discrete ordinate method described in Chapter 11 of [2]. In particular, we use the Henyey–Greenstein scattering phase function to compute the redistribution function defined as 21 − g2 2b p E ; (6.2) ab πa − b a b pp where a 1 g2 − 2gμμ0 , b 2g 1 − μ2 1 − μ02 , and Ek is the complete elliptic integral of the second kind. Upon numerical solution of boundary value problem (6.1), we compute F defined in Eq. (4.1) using the same Gauss-Legendre quadrature rule used for the discrete ordinate method. For the gKM system, we calculated numerically the entries of H 1 and H 2 using h defined in Eq. (6.2) and GaussLegendre quadrature. We also compute the vectors f defined in Eqs. (3.10) and (3.11) using Gauss-Legendre quadrature. With those matrices and vectors defined, we solve numerically the generalized eigenvalue problem (4.7) and obtain λ~ j , u~ j , and v~ j for j 1; 2. The solution of the gKM equations ~ is given by Eq. (4.9) with a~ j τ eλj τ−τ2 α~ j , and b~ j τ hμ; μ0
e−λj τ−τ1 β~ j . Evaluating this solution at τ τ1 and substituting that result into Eq. (4.6a), we find that ~
2 n X
o ~ u~ j e−λj τ2 −τ1 α~ j v~ j β~ j Mf :
(6.3)
j1
Evaluating the solution at τ τ2 and substituting that result into Eq. (4.6b), we find that
C. Sandoval and A. D. Kim 2 n X
Vol. 31, No. 3 / March 2014 / J. Opt. Soc. Am. A
o ~ v~ j α~ j u~ j e−λj τ2 −τ1 β~ j Mf − :
(6.4)
j1
Combining Eqs. (6.3) and (6.4) yields a 4 × 4 linear system of equations for α~ j and β~ j for j 1; 2. With α~ j and β~ j determined, we compute F through evaluation of F τ
2 n o X ~ ~ U~ 1j eλj τ−τ2 α~ j V~ 1j e−λj τ−τ1 β~ j ;
(6.5)
j1
F − τ
2 n X
o ~ ~ V~ 1j eλj τ−τ2 α~ j U~ 1j e−λj τ−τ1 β~ j :
(6.6)
j1
Here, U~ 1j denotes the first entry of u~ j , and V~ 1j denotes the first entry of v~ j . We solve the KM system (1.1) using K and S defined in Eqs. (5.10) and (5.11), respectively. We supplement this system with the boundary conditions Z1 F τ1 f μμdμ; (6.7) 0
F − τ2
Z
1
0
f − −μμdμ:
(6.8)
633
We compute the integrals in Eqs. (6.7) and (6.8) using GaussLegendre quadrature. In Fig. 1 we show results comparing solutions to (6.1) computed using the RTE (solid curves), the gKM (dashed curves), and the KM equations (dot-dashed curves). Here, we have set τ1 0, τ2 1, ϖ 0 0.99, and g 0. The top plot shows F as a function of τ, and the bottom plot shows F − as a function of τ. The gKM theory gives a more accurate solution compared to KM theory. It is nearly indistinguishable from that computed through the solution of the RTE. In Fig. 2 we plot the results when ϖ 0 0.99 and g 0.8. For this case, scattering is forward-peaked. Consequently, the specific intensity is not approximated well by a piecewise linear approximation. Thus, the DP 1 is not as good of an approximation, which, in turn, means that the gM and KM equations are less accurate. This is evidenced by the fact that both the gKM and KM results are less accurate than those shown in Fig. 1. Nonetheless, Fig. 2 shows that the gKM equations are still much more accurate than the KM equations. Next, we investigate reducing the value of ϖ 0 . For this case, the KM equations are known to be inaccurate [8,16,18]. Figure 3 shows results for ϖ 0 0.5 and g 0. Although scattering is isotropic, so that g 0, we find that the KM equations are much less accurate than the gKM equations. The error is worse for F τ than for F − τ presumably because it must also take into account the nonconstant boundary condition
0.3 RTE gKM KM
0.25
0.18 RTE gKM KM
0.16
F+(τ)
0.2
+
F (τ)
0.14
0.15
0.12
0.1 0.1
0.05
0
0.2
0.4
0.6
0.8
1
0.08
τ
0
0.2
0.4
0.55 RTE gKM KM
0.5
τ
0.6
0.8
1
0.6
0.8
1
0.52 RTE gKM KM
0.5
0.45
−
F (τ)
−
F (τ)
0.48
0.4
0.46 0.44
0.35 0.42 0
0.2
0.4
0.6
0.8
1
τ
F τ
F − τ
Fig. 1. Comparisons of (top) and (bottom) computed from solutions of boundary value problem (6.1) using the radiative transport equation (RTE), the generalized Kubelka–Munk (gKM) equations, and the Kubelka–Munk (KM) equations. Here, τ1 0, τ2 1, ϖ 0 0.99, and g 0.
0.4
0
0.2
0.4
τ Fig. 2. Comparisons of F τ (top) and F − τ (bottom) computed from solutions of boundary value problem (6.1) using the RTE, the gKM equations, and the KM equations. Here, τ1 0, τ2 1, ϖ 0 0.99, and g 0.8.
634
J. Opt. Soc. Am. A / Vol. 31, No. 3 / March 2014
C. Sandoval and A. D. Kim
(6.1b). When scattering is also forward-peaked, as shown in Fig. 4, we find that the KM equations produce negative values for F , which are unphysical. However, even with ϖ 0 0.5, the gKM equations yield accurate results. Even at a lower albedo (ϖ ∼ 0), the gKM equations yield qualitatively correct results. Although we do not show these results here, we have found that computing solutions with smaller (τ2 − τ1 ∼ 0.1) and larger (τ2 − τ1 ∼ 10) optical thicknesses yield similar results to those shown here. B. Nonhomogeneous Problem for the Diffuse Flux We now study the boundary value problem μ
∂I ϖ I 0 ∂τ 2
Z
1
−1
hμ; μ0 Iμ0 ; τdμ0
in τ1 < τ < τ2 ;
(6.9a)
Iμ; τ1 δμ − 1 on 0 < μ ≤ 1;
(6.9b)
Iμ; τ2 0 on − 1 ≤ μ < 0:
(6.9c)
This boundary value problem governs a plane wave with unit incident flux normally on a plane-parallel slab. It is a fundamental problem in the study of radiative transfer [1,2]. A complication in solving boundary value problem (6.9) is the Dirac delta function appearing in boundary condition (6.9b). To study this problem, one represents the specific intensity as
the sum I I ri I d . Here, I ri denotes the reduced intensity satisfying μ
(6.10a)
I ri μ; τ1 δμ − 1 on 0 < μ ≤ 1;
(6.10b)
I ri μ; τ2 0 on − 1 ≤ μ < 0:
(6.10c)
The solution of boundary value problem (6.10) is readily found to be I ri δμ − 1e−τ−τ1 ∕μ . The diffuse intensity, I d , satisfies the following nonhomogeneous, boundary value problem for the RTE: μ
∂I d ϖ Id 0 ∂τ 2
Z
1
−1
hμ; μ0 I d μ0 ; τdμ0 Qri ;
(6.11a)
I d μ; τ1 0 on 0 < μ ≤ 1;
(6.11b)
I d μ; τ2 0 on − 1 ≤ μ < 0:
(6.11c)
The nonhomogeneous term, Qri 1∕2ϖ 0 hμ; 1e−τ−τ1 , in Eq. (6.11a) comes from inserting I ri into the integral operation in the RTE. 0.1
0.11 RTE gKM KM
0.1
0.08 0.06 0.04
+
+
F (τ)
F (τ)
0.09
0.02
0.08
0 −0.02
RTE gKM KM
0.07 −0.04 0.06
∂I ri I ri 0; ∂τ
−0.06 0
0.2
0.4
τ
0.6
0.8
0
1
0.2
0.4
τ
0.6
0.8
1
0.6
0.8
1
0.55
0.55 RTE gKM KM
0.5
RTE gKM KM
0.5
0.45
0.45 F (τ)
0.35
0.4
−
−
F (τ)
0.4
0.35
0.3
0.3
0.25
0.25
0.2
0.2 0
0.2
0.4
0.6
0.8
1
τ Fig. 3. Comparisons of F τ (top) and F − τ (bottom) computed from solutions of boundary value problem (6.1) using the RTE, the gKM equations, and the KM equations. Here, τ1 0, τ2 1, ϖ 0 0.5, and g 0.
0
0.2
0.4 τ
Fig. 4. Comparisons of F τ (top) and F − τ (bottom) computed from solutions of boundary value problem (6.1) using the RTE, the gKM equations, and the KM equations. Here, τ1 0, τ2 1, ϖ 0 0.5, and g 0.8.
C. Sandoval and A. D. Kim
Vol. 31, No. 3 / March 2014 / J. Opt. Soc. Am. A
results for y y H yP at τ τ1 and substituting that result into Eq. (4.6a) with f 0, we find that
We compute the diffuse fluxes, defined as F d τ
Z
1 0
I d μ; τμdμ;
(6.12)
by solving boundary value problem (6.11) directly and by solving the gKM for this problem. To solve boundary value problem (6.11), we use the same method that we use to solve boundary value problem (6.1). The gKM equations for this problem take the form
I 0 d y S S 1 2 0 −I dτ y− S2 S1
ϖ 0 h −τ−τ1 y e ; (6.13) y− 2 h−
where the entries of the vectors h are given by Z1 h hμ; 1P n μdμ; n 0; 1: n1
−I − S 1 −S 2
−S 2 I − S1
η η−
ϖ h 0 − : 2 h
2 n X
o ~ v~ j α~ j u~ j e−λj τ2 −τ1 β~ j −η− e−τ2 −τ1 :
(6.17)
j1
Combining Eqs. (6.3) and (6.4) yields a 4 × 4 linear system of equations for α~ j and β~ j for j 1; 2. With α~ j and β~ j determined, we compute F through evaluation of 2 n X
o ~ ~ −τ−τ1 ; (6.18) U~ 1j eλj τ−τ2 α~ j V~ 1j e−λj τ−τ1 β~ j η 1 e
j1
F − τ
2 n o X ~ ~ V~ 1j eλj τ−τ2 α~ j U~ 1j e−λj τ−τ1 β~ j η−1 e−τ−τ1 ;
(6.19) with η 1 denoting the first entries of η . In Fig. 5, we compare the diffuse fluxes computed from solutions of (6.11) by solving the RTE (solid curves) and
0.05
0.1
0.04
0.08
0.03
0.06
F (τ)
0.02
+
+
Evaluating the solution at τ τ2 and substituting that result into Eq. (4.6b) with f − 0, we find that
(6.15)
The homogeneous solution is given by Eq. (4.9) with ~ ~ a~ j τ eλj τ−τ2 α~ j , and b~ j τ e−λj τ−τ1 β~ j . Evaluating our
0.01
0.04 0.02
0
0
RTE gKM 0
0.2
0.4
τ
0.6
0.8
−0.02
1
0.06
0.016
0.05
0.014
RTE gKM 0
0.2
0.4
τ
0.6
0.8
1
0.6
0.8
1
0.012
0.04
0.01
F (τ)
0.03
−
−
F (τ)
(6.16)
j1
y H
−0.01
o ~ u~ j e−λj τ2 −τ1 α~ j v~ j β~ j −η :
j1
F τ
We compute the solution of Eq. (6.13) as the sum of the homo geneous and particular solutions: y y H yP . We seek the −τ−τ 1 , where η particular solution in the form yP η e satisfies the 4 × 4 linear system
F (τ)
2 n X
(6.14)
0
635
0.02
0.008 0.006
0.01
0.004
0 −0.01
RTE gKM 0
0.2
0.002 0.4
0.6
0.8
1
τ Fig. 5. Comparisons of F τ (top) and F − τ (bottom) computed from solutions of boundary value problem (6.11) using the RTE, and the gKM equations. Here, τ1 0, τ2 1, ϖ 0 0.99, and g 0.
0
RTE gKM 0
0.2
0.4
τ
Fig. 6. Comparisons of F τ (top) and F − τ (bottom) computed from solutions of boundary value problem (6.11) using the RTE and the gKM equations. Here, τ1 0, τ2 1, ϖ 0 0.99, and g 0.7.
636
J. Opt. Soc. Am. A / Vol. 31, No. 3 / March 2014
the gKM equations (dashed curves). Here, we have set τ1 0, τ2 1, ϖ 0 0.99, and g 0. Here, the gKM equations provide a very accurate approximation to the RTE. Figure 6 compares the diffuse fluxes for the case in which g 0.7. As expected, the approximation produced by the gKM equations are less accurate when scattering is forward-peaked than when it is isotropic. Nonetheless, the gKM equations still provide an accurate approximation to the RTE for the diffuse fluxes. Although we do not show these results here, we find that gKM equations still approximate the RTE in the manner shown in Figs. 5 and 6 with lower albedo values, and for both smaller and larger optical thicknesses.
7. CONCLUSION We have systematically derived the KM equations from the RTE. This derivation involves the application of the DP N method to solve the RTE. By transforming the DP 1 system of equations, we obtained the gKM equations. Then, by studying the DP 1 system of equations in the limit of strong multiple scattering, we derived the KM equations. For the KM equations, we have determined all coefficients explicitly. In particular, we have found that K~ 2μa and S~ 3∕4μs 1 − g − 1∕4μa 1–3g. These results agree closely with those found in the literature. However, the significant difference in the results found here is the term in S~ proportional to μa g. Moreover, we established this theory’s range of validity as a consequence of the asymptotic analysis used to derive it. The gKM equations can easily account for general boundary conditions and nonhomogeneous problems. The numerical results shown here demonstrate that the gKM equations are much better than the KM equations at approximating the solutions of the RTE. The complexity of the gKM equations is only slightly larger, since it is a 4 × 4 system rather than a 2 × 2 one. Since it is more accurate than the KM system of equations, more broadly applicable, and only slightly more difficult to solve, the gKM equations should be very useful.
REFERENCES 1. 2.
S. Chandrasekhar, Radiative Transfer (Dover, 1960). A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE, 1996).
C. Sandoval and A. D. Kim 3. P. Kubelka and F. Munk, “Ein beitrag zur optik der farbanstriche,” Z. Tech. Phys. 12, 593–601 (1931). 4. P. Kubelka, “New contributions to the optics of intensely lightscattering materials. Part I,” J. Opt. Soc. Am. 38, 448–457 (1948). 5. B. Philips-Invernizzi, D. Dupont, and C. Cazé, “Bibliographical review for reflectance of diffusing media,” Opt. Eng. 40, 1082–1092 (2001). 6. B. J. Brinkworth, “Interpretations of the Kubelka-Munk coefficients in reflection theory,” Appl. Opt. 11, 1434–1435 (1972). 7. L. F. Gate, “Comparison of the photon diffusion model and Kubelka-Munk equation with exact solution of the radiative transport equation,” Appl. Opt. 13, 236–238 (1974). 8. J. Nobbs, “Kubelka-Munk theory and the prediction of reflectance,” Rev. Prog. Coloration 15, 66–75 (1985). 9. W. M. Star, J. P. A. Marijnissen, and M. J. C. Van Gemert, “Light dosimetry in optical phantoms and in tissues: I. Multiple flux and transport theory,” Phys. Med. Biol. 33, 437–454 (1988). 10. W. E. Vargas and G. A. Niklasson, “Applicability conditions of the Kubelka-Munk theory,” Appl. Opt. 36, 5580–5586 (1997). 11. R. Molenaar, J. ten Bosch, and J. Zijp, “Determination of Kubelka-Munk scattering and absorption coefficients by diffuse illumination,” Appl. Opt. 38, 2068–2077 (1999). 12. L. Yang and B. Kruse, “Revised Kubelka-Munk theory. I. Theory and application,” J. Opt. Soc. Am. A 21, 1933–1941 (2004). 13. L. Yang, B. Kruse, and S. J. Miklavcic, “Revised Kubelka-Munk theory. II. Unified framework for homogeneous and inhomogeneous optical media,” J. Opt. Soc. Am. A 21, 1942–1952 (2004). 14. L. Yang and S. J. Miklavcic, “Revised Kubelka-Munk theory. III. A general theory of light propagation in scattering and absorptive media,” J. Opt. Soc. Am. A 22, 1866–1873 (2005). 15. P. Edström, “Examination of the revised Kubelka-Munk theory: considerations of modeling strategies,” J. Opt. Soc. Am. A 24, 548–556 (2007). 16. S. N. Thennadil, “Relationships between the Kubelka-Munk scattering and radiative transfer coefficients,” J. Opt. Soc. Am. A 25, 1480–1485 (2008). 17. M. Neuman and P. Edström, “Anisotropic reflectance from turbid media. I. Theory,” J. Opt. Soc. Am. A 27, 1032–1039 (2010). 18. M. L. Myrick, M. N. Simcock, M. Baranowski, H. Brooke, S. L. Morgan, and J. N. McCutcheon, “The Kubelka-Munk diffuse reflectance formula revisited,” Appl. Spectrosc. Rev. 46, 140–165 (2011). 19. B. Davison, Neutron Transport Theory (Oxford University, 1958). 20. K. M. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, 1967). 21. E. E. Lewis and W. F. Miller, Computational Methods of Neutron Transport (American Nuclear Society, 1993). 22. R. Aronson, “PN vs. double-PN approximations for highly anisotropic scattering,” Transp. Theory Stat. Phys. 15, 829–840 (1986).
