A Fast Algorithm to the Radiative
Transport Equation and Implementationof Theory Into an Applet
Mihai V. Avram(Purdue University)
Faculty Mentor: Manabu Machida(University of Michigan)
7/24/2013

Abstract
In random media, such as clouds or biological tissue, light obeys theRadiative Transport Equation (RTE). Using the RTE, we developed a fastalgorithm and implemented a Java applet to calculate the energy density oflight in three dimensions.

1

Introduction

The radiative transport equation (RTE) has a connection to many interdisciplinary fields [1, 2, 3, 4]. In particular, light in random media is
described by the RTE. For example, optical tomography is formulated as
inverse problems of the RTE [5, 6].
A lot of numerical methods of solving the RTE have been developed[1, 2, 3, 4, 5, 6, 7], which include the Monte Carlo method, the finite element
method, the PL method, and the method of discrete ordinates. In 1960,Professor Kenneth Myron Case from University of Michigan published apaper on solving the one-dimensional RTE analytically [8]. Our theory inthis report for the three-dimensional RTE is based on Case’s method [9].
In this paper we will describe the theory for the three dimensional radiative transport equation. We will also delve into the implementation of aJava applet. The applet is programmed to graph the energy density of lightin a three dimensional random medium.
1

2

Solving the Radiative Transport Equation

We consider the time-independent radiative transport equation (RTE)
in three dimensions:Zp(ŝ, ŝ0 )u(~r, ŝ)dŝ0 + S(~r, ŝ).
(1)ŝ · ∇u + (µa + µs )u = µs
S2Here variables are defined as follows:~r: The three dimensional position vector ~r = hx, y, zi.
ŝ: The three dimensional directional vector ŝ = hsin(θ) cos(φ), sin(θ)sin(φ), cos(θ)i.
∂∂∂∇: The three dimensional gradient vector ∇ = h ∂x, ∂y, ∂zi.
u(~r, ŝ): The intensity of light dependent on position ~r and direction ŝ.
µa : The absorption constant proportional to the probability ofabsorption per unit length.
µs : The scattering constant proportional to the probability of scattering
per unit length.
S(~r, ŝ): Source of light.
p(ŝ, ŝ0 ): The phase scattering function, modeling the probability of lightbeing scattered from direction ŝ0 to direction ŝ.
Here we assume linear scattering
13g+(ŝ · ŝ0 ), g ∈ [0, 1].
4π 4πNote that g is the linear asymmetry parameterZg=(ŝ · ŝ0 )p(ŝ, ŝ0 )dŝ0 .
p(ŝ, ŝ0 ) =(2)

(3)

S2

2.1

Specific Intensity
We define the total attenuation µt asµt = µs + µa .

(4)

Let us now introduce a constant c:
c=µs.
µa + µ s(5)By taking the unit of length to be 1/µt , we divide (1) by µt . Thus, werewrite (1) into a more convenient form:
ZS(~r, ŝ)ŝ · ∇u + u = cp(ŝ, ŝ0 )u(~r, ŝ)dŝ0 +.
(6)
µt
S22

We then obtain u(~r, ŝ) as a superposition of elementary solutions uν (~r, ŝ; ~q),which obey
Zp(ŝ, ŝ0 )uν (~r, ŝ0 )dŝ0 .

ŝ · ∇uν + uν = c(7)S2The elementary solutions are labeled by ν ∈ R and ~q ∈ R2 and are obtainedby separation of variables [9].

2.2Eigenvalues

In the case of linear scattering (2), the index ν which labels uν , is either±ν0 where (ν0 > 1), or any value on (−1, 1). We can calculate ν0 as follows[9, 10]. Let us expand uν with Legendre polynomials.
∞ rX
2l + 1uν (~r, ŝ; ~q) =Cl Pl (ŝ · k̂)e−k̂·~r/ν ,(8)4πl=0where the vector k̂ is given by [9]p−iν~q, Q(νq) = 1 + ν 2 q 2 ,k̂ =Q(νq)

q = |~q| .

(9)Here Legendre polynomials are recursively computed as
P0 = 1,P1 = µ,Pl+1 =(2l + 1)µPl − lPl−1.
l+1

Moreover, spherical harmonics are given bys
2l + 1 (l − m)! mYlm (θ, ϕ) =P (cos θ)eimφ .
4π (l + m)! l(10)(11)

By plugging (8) into (7), we end up with the following equation.

ŝ · ∇∞
X
l=0r

2l + 1Cl Pl (ŝ · k̂)e−k̂·~r/ν
4π

∞X

r

2l + 1Cl Pl (ŝ · k̂)e−k̂·~r/ν
4πl=0Z∞ r
X
2l + 10= cp(ŝ, ŝ )Cl Pl (ŝ0 · k̂)e−k̂·~r/ν dŝ0 .
4πS2+l=0(12)3We writep(ŝ, ŝ0 ) =X∗pl Ylm (ŝ · k̂)Ylm(ŝ0 · k̂),(13)lmwhere p0 = 1, p1 = g, p2 = p3 = · · · = 0. We now haver∞ r∞ rXXµ X 2l + 12l + 12l + 1−Cl Pl (µ) +Cl Pl (µ) = cpl
Cl Pl (µ).
ν4π
4π4πl=0ll=0(14)We then multiply Pl0 (µ) and integrate both sides over µ. Note thatZ 1Pl (µ)Pl0 (µ)dµ =−12δll0 ,2l + 1(15)where δll0 is the Kronecker’s delta function (if l = l0 , then δll0 = 1 and if
l 6= l0 , then δll0 = 0). After some calculations, (14) reduces to the following
equation:
#"∞X
√
√
l+1lp
δl0 ,l+1 + p
δl0 ,l−1 Cl 2l + 1 σl2
2[4(l + 1) − 1](σl+1 )(σl )(4l − 1)(σl−1 )(σl )l=0p√= νCl0 2l0 + 1 σl0 ,(16)where σl = (1 − cpl ).
Upon rewriting (14) in the form of a matrix-vector equation Bx = νx,
we can solve for eigenvalues ν and ν0 as explained below.

B00 B01 B02 .
.
B0n
 B10 B11 B12 .
.
B1n 

 B20 B21 B22 .
.
B2n
.
B=
 .
.
.
.
.
. 


 .
.
.
. Bn−1n−1. 
Bn0 Bn1 Bn3 .

.

Bnn

Note that from (16) we have#"ll+1
δl0 ,l+1 + p
δl0 ,l−1 .
Bl0 l = p[4(l + 1)2 − 1](σl+1 )(σl )(4l2 − 1)(σl−1 )(σl )(17)4That is,0 √ 1
 3σ1 σ0

0

B=0

0

0

√ 13σ1 σ00000.
000.
0.

000

√ 215σ2 σ10√ 215σ2 σ100.
.
0√n(4n2 −1)σn σn−10000






.


√ 2 n

(4n −1)σn σn−1 0Hence, diagonalizing the symmetric matrix B will yield real eigenvaluesν ∈ (0, 1) and ν0 > 1, where ν0 is the largest eigenvalue. These eigenvalues
can later be used to solve for constants N0 and N (ν).

3Finding Energy Density U (z)The Green’s function G(~r, ŝ; r~0 , ŝ0 ) for the radiative transport equation
satisfies
Zŝ · ∇G + G = cp(ŝ, ŝ0 )G(~r, ŝ0 )dŝ0 + δ(~r − r~0 )δ(ŝ − ŝ0 ),(18)
S2

where δ(~r − r~0 ), δ(ŝ − ŝ0 ) are the Dirac delta functions with properties:
Z ∞δ(x)dx = 1,(19)−∞Z ∞f (x)δ(x − x0 )dx = f (x0 ).

(20)

−∞

We obtain G in (18) using uν in (7), which we write as uν (~r, ŝ; ~q) = φν (µ(k̂))e−k̂·~r/ν ,where µ(k̂) = ŝ · k̂. We can write G asG+ = A+ φν0 (µ(k̂))eZ 1i~q ·~ρ −Q(ν0 q)z/ν0e+A(ν)φν (µ(k̂))ei~q·~ρ e−Q(νq)z/ν dν

0(21)for z > z0 ,i~q ·~ρ Q(ν0 q)z/ν0

G− = −A− φ−ν0 (µ(k̂))ee

Z 0−A(ν)φν (µ(k̂))ei~q·~ρ e−Q(νq)z/ν dν−1(22)
5for z < z0 , where ρ~ = hx, yi. Moreover, A+ , A− , and A(ν) are functions of~q that we must solve Rfor. To solve for A+ , A− , and A(ν), we integrate (18)z +with respect to z as z00− (18)dz for  ≈ 0. Then we use the orthogonality
relations for φν (µ):
Z 1

µφ2ν0 (µ)dµ = N0 ,(23)µφν (µ)φ0ν (µ)dµ = N (ν)δ(ν − ν 0 ),(24)−1

Z 1−1Z 1−1µφν (µ)φ0ν (µ)dµ = 0,(25)if ν 6= ν 0 . Here N0 and N (ν) are given by"1(3γ(ν0 ) − 2)(1 − c + cγ(ν0 ))cν03c
ν2N0 =γ(1) − 0γ(ν0 ) 22γ(ν0 )ν0 − 1#+ 3c(γ(1) − 1) ,(26)ihπcνN (ν) = ν (1 − c + cγ(ν) − cνγ(ν) tanh−1 (ν))2 + (γ(ν))2 ,2(27)where γ(ν) = 1 + 3gσ0 ν 2 . Using the orthogonality relations and some Diracdelta function relations we obtain:
1e−iq~0 ·ρ~0 eQ(ν0 q0 )z0 /ν0 φ∗ν0 (µ0 (k̂)),
2πQ(ν0 q0 )N0−1
A− (q~0 ) =e−iq~0 ·ρ~0 e−Q(ν0 q0 )z0 /ν0 φ∗−ν0 (µ0 (k̂)),
2πQ(ν0 q0 )N01e−iq~0 ·ρ~0 eQ(νq0 )z0 /ν φ∗ν (µ0 (k̂)).
A(ν) =2πQ(νq0 )N (ν)
A+ (q~0 ) =(28)
(29)
(30)

We can then include A+ , A− , and A(ν) in (21) and (22) and combining (21)and (22) will yield:"Z∗−Q(ν0 q)|z−z0 |1i~q ·(~ρ−ρ~0 ) φ±ν0 (µ(k̂))φ±ν0 (µ0 (k̂))νG(~ρ, z, ŝ; ρ~0 , z0 , ŝ0 ) =ee(2π)3 R2Q(ν0 q)N0#Z 1
φ±ν0 (µ(k̂))φ∗±ν0 (µ0 (k̂)) −Q(νq)|z−z0 |ν+edν d~q.
(31)
Q(νq)N (ν)
06We now recover the dimension. Let us consider a point sourceS(~r, ŝ) = δ(~r).

(32)

We obtain u asu(~ρ, z, ŝ) = µ2t

ZZG(~ρ, z, ŝ; ρ~0 , z0 , ŝ0 )δ(ρ~0 )δ(z0 )d~ρ0 dz0 dŝ0 .

R3(33)

S2

We calculate the energy density U of light defined asZ
1
U=u(~r, ŝ)dŝ,v(34)

where v is the speed of light in the medium. Since U is spherically symmetric,
we measure it along the z-axis as −µ z
Z 1 ( −µt z )( νt )µt  e 0e νU (z) =dν  , z > 0.
(35)+vz ν0 N00 νN (ν)
We conclude the theory for our algorithm by plotting U (z) in Figure 1.

7

Figure 1: Energy density U (z) in (35) for different optical parameters compared with results from Monte Carlo simulations [13]. For isotropic scattering, U (z) can also be obtained analytically with the Fourier transform.

4

Implementation of the Applet

We aspired to create an applet, in the programming language Java,
which can graph the intensity of light depending on a point in three dimensional space.
A snippet of our applet can be seen below:

8Figure 2: Java Applet Using Quick RTE Algorithm.
The applet simply takes the user input: µa , µs , and g (see page 2 ofthis report) which are the absorption constant, scattering constant, andscattering asymmetry parameter, respectively. Upon invoking the ”Graph”
button, a list of [z, U (z)] coordinates appear on the list to the right of thegraphic, and a density line is drawn on the graph. The coordinates in the1list are exact measures of intensity accurate to 1000
th of a decimal.
In the graph, the vertical axis shows (ln U )sc , which is given below, andthe horizontal axis is z/lt .
We discretize z as
zi = i∆z,

i = 1, 2, ..., N,

N = 50,N ∆z = 10lt ,(36)where lt = 1/µt . We then define (ln U )sc as(ln U )sc =ln U (zi ) − ln U (zN ),ln U (z1 ) − ln U (zN )(37)This scaling prevents the graph from digressing from the viewing window ofz/lt ∈ [0, 10] and (ln U )sc ∈ [0, 1] in our applet.

95Summary

In this report we illustrated the theory for solving the three dimensional radiative transport equation. Moreover, we implemented a Java applet which is formulated to graph the energy density of light in a three
dimensional random medium. The URL for the Java applet can be foundbelow.
URL: https://sites.google.com/site/ezrte13/

References[1] S. Chandrasekhar, Radiative Transfer (Dover, 1960).
[2] K. M. Case and P. F. Zweifel, Linear Transport Theory (AddisonWesley, 1967).
[3] A. Ishimaru Wave Propagation and Scattering in Random Media (Academic Press, 1978).
[4] L. A. Apresyan and Y. A. Kravtsov, Radiation Transfer: Statistical andWave Aspects (Gordon and Breach Publishers, 1996).
[5] S. R. Arridge “Optical tomography in medical imaging,” Inverse Problems 15 (1999) R41–R93.
[6] S. R. Arridge and J.C. Schotland, “Optical Tomography: Forward andInverse Problems,” Inverse Problems 25 (2009) 123010 (pp. 59).
[7] K. Ren, “Recent Developments in Numerical Techniques for TransportBased Medical Imaging Methods,” Commun. Comput. Phys. 8 (2010)
1–50.
[8] K. M. Case “Elementary Solutions of the Transport Equation and TheirApplications,” Ann. Phys. 9 (1960) 1–23.
[9] M. Machida, “The three-dimensional radiative transport equation with
Case’s method and rotated reference frames,” arXiv:1304.1941.
[10] M. Machida,
arXiv:1304.5207.

“Case’sEigenvalues

byMarkel’sMatrix,”[11] V. A. Markel, “Modified spherical harmonics method for solving theradiative transport equation,” Waves in Random Media 14 (2004) L13–L19.
10

[12] G. Panasyuk, J. C. Schotland, and V. A. Markel, “Radiative transport
equation in rotated reference frames,” J. Phys. A: Math. Gen. 39 (2006)
115–137.
[13] The Monte Carlo package MC written by Vadim A. Markel was used
after implementing linear scattering.

11