Research Article

Vol. 33, No. 4 / April 2016 / Journal of the Optical Society of America A

551

Radiative transport and optical tomography withlarge datasets
MANABU MACHIDA,1 GEORGE Y. PANASYUK,2 ZHENG-MIN WANG,3 VADIM A. MARKEL,4 AND
JOHN C. SCHOTLAND1,5,*
1Department of Mathematics, University of Michigan, Ann Arbor, Michigan 48109, USA
Propulsion Directorate, Air Force Research Laboratory, Wright-Patterson Air Force Base, Ohio 45433, USA
3
Department of Bioengineering, University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA
4
Aix-Marseille Université, CNRS, Centrale Marseille, Institut Fresnel UMR 7249, 130113 Marseille, France
5
Department of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA*Corresponding author: schotland@umich.edu
2Received 21 September 2015; revised 28 October 2015; accepted 16 November 2015; posted 11 December 2015 (Doc. ID 250385);
published 10 March 2016

We consider the inverse problem of optical tomography in the radiative transport regime. We report numericaltests of a direct reconstruction method that is suitable for use with large datasets. Reconstructions of experimental
data obtained from a noncontact optical tomography system are also reported. © 2016 Optical Society of America
OCIS codes: (170.6960) Tomography; (110.0113) Imaging through turbid media.
http://dx.doi.org/10.1364/JOSAA.33.000551

1. INTRODUCTION
The development of tools to probe the structure of highly scattering media such as clouds, colloids, and biological tissue is offundamental interest and considerable applied importance.
One such method, known as optical tomography (OT), is abiomedical imaging modality with unique capabilities to assess
physiological function, such as blood volume and tissue oxygenation [1]. In a typical OT experiment, a material medium
is illuminated by a narrow collimated beam and the light that
propagates through the medium is collected by an array ofdetectors. In many instruments, the sources and detectors
are coupled to the medium by means of optical fibers andthe number of measurements which can be obtained varies
from 102 − 105 source-detector pairs. More recently, noncontact
imaging systems have been introduced, wherein a scanned
beam and a lens-coupled CCD camera are employed toreplace the aforementioned illumination and detection fiberoptics [2–4]. Using such a noncontact method, extremely large
datasets of order 107 − 109 measurements can readily be
acquired. The availability of such datasets has been shownto vastly improve the quality of images in OT [4–6].
The inverse problem of OT is to reconstruct the opticalproperties of a highly scattering medium from boundary
measurements. The standard computational approach to thisproblem is to minimize a penalized least-squares functional[1,7–15]. Although such optimization methods are extremely
flexible, they have very high computational cost and are
not well suited to the large datasets of noncontact optical
1084-7529/16/040551-08 Journal © 2016 Optical Society of America

tomography. Direct reconstruction methods offer an alternative
approach to optimization-based algorithms [16–23]. By directreconstruction, we mean the use of inversion formulas andassociated fast algorithms. Such formulas have been derived
for particular experimental geometries, including those with planar boundaries. To date, these methods have been developed
within the framework of the diffusion approximation (DA)to the radiative transport equation (RTE). We note that the
DA is accurate when the energy density of the optical field varies
slowly on the scale of the transport mean free path. The DA
breaks down in optically thin layers, near boundaries and instrongly absorbing media, conditions that are frequently
encountered in biomedical applications.
In this paper, we develop a direct reconstruction method forthe inverse problem in the transport regime. The primary toolsthat we exploit are an integral equation formulation of the forward problem [24] and the recently derived plane-wave decomposition for the Green’s function of the RTE [25–27]. Using
this approach, we show that the linearized inverse problem canbe formulated in terms of the inversion of a suitably defined
Fourier–Laplace transform. We illustrate this result withnumerical simulations and reconstructions of experimental dataobtained from a noncontact OT system.
The remainder of this paper is organized as follows.
In Section 2, we introduce the linearized forward problem
within the accuracy of the Rytov approximation. This is followed in Section 3 by the derivation of the Fourier–Laplace
inversion formula for the slab geometry. Numerically simulated

552

Vol. 33, No. 4 / April 2016 / Journal of the Optical Society of America A

reconstructions are presented in Section 4. Reconstructions
from experimental data are shown in Section 5. Our conclusions are presented in Section 6. The derivation of theGreen’s function for the RTE with reflecting boundary conditions is given in Appendix A.
2. FORWARD PROBLEM

Research Article

It will prove convenient to introduce the data function ϕ,which is defined byϕ  −I 0 lnI ∕I 0 :

(9)It follows from Eq. (8) that ϕ obeys the integral equationZ(10)
ϕr  d3 r 0 d2 s 0 Gr; ŝ; r 0 ; ŝ 0 δμa r 0 I 0 r 0 ; ŝ 0 :

A. Radiative Transport

The propagation of multiply scattered light in a volume Ω istaken to be governed by the RTE:
Z
ŝ·∇I μa rμs rI  μs r Aŝ; ŝ 0 I r; ŝ 0 d2 s 0 Sr;ŝ:(1)Here, I r; ŝ denotes the specific intensity at the point r in the
direction ŝ, S is the source, and μa and μs are the absorption andscattering coefficientsR of the medium. The phase function A isnormalized so that Aŝ; ŝ 0 d2 s 0  1 for all ŝ. The specificintensity is also assumed to obey the half-range boundary
conditionI r; ŝ 0   Rjŝ · ŝ 0 jI r; ŝ;n̂ · ŝ 0 < 0 on ∂Ω:(2)Here, R is the reflection coefficient, n̂ is the outward unit normal to ∂Ω, and ŝ and ŝ 0 are the incident and reflected
directions at the boundary, respectively. See Appendix A for
further details.
We assume that the scattering coefficient μs is constanteverywhere in the medium but the absorption coefficient μa
varies with position. We thus decompose μa into a constantpart μ̄a and a spatially varying part δμa :μa r  μ̄a  δμa r:The RTE (1) thus becomesŝ · ∇I  μt I  δμa rI  μs

Z

0

(3)

0

2 0

Aŝ; ŝ I r; ŝ d s  S;

(4)

where μt  μ̄a  μs . The solution to Eq. (4) is given byZ
I r; ŝ  I 0 r; ŝ − d3 r 0 d2 s 0 Gr; ŝ; r 0 ; ŝ 0 δμa r 0 I r 0 ; ŝ 0 :
(5)Here, I 0 is the incident-specific intensity, defined as thesolution to Eq. (1) with μa  μ̄a obeying the boundary condition (2). The Green’s function G satisfies the equationŝ · ∇r Gr; ŝ; r 0 ; ŝ 0   μt Gr; ŝ; r 0 ; ŝ 0  Zμs

B. Slab Geometry

For the remainder of this paper we consider a threedimensional slab-shaped medium, as shown in Fig. 1. The incident field is taken to be generated by a point source oriented
in the inward normal direction located on the z  0 plane.
Thus I 0  Gr; ŝ; ρs ; 0; ẑ, where ρs is the transverse coordinate of the source. Light exiting the slab in the outward normal
direction is collected by a point detector that is located on theplane z  L. The integral equation (10) thus becomesZ
ϕρs ; ρd   d3 rd2 sGρd ; L; ẑ; r; ŝGr; ŝ; ρs ; 0; −ẑδμa r;
(11)
where ρd is the transverse coordinate of the detector.
To make further progress, we make use of the method ofrotated reference frames to express the Green’s function as
an expansion in two-dimensional plane waves and sphericalharmonics [25–27]. The details are presented in Appendix A,
where it is shown thatGr; ŝ;r 0 ; ŝ 0 
Z 2
∞ X
l
d q −iq·ρ−ρ 0  Xk l m q;zY l m ŝY l m ŝ 0 ;2e2πl 0 m−lwhere r  ρ; z and kl m is defined in Eq. (A19). Using thisresult, we find that Eq. (11) becomesZ
Z 2 Z 2
d q1
d q2 iq −q ·ρ −iq ·ρs −q ·ρd 
3ϕρs ; ρd   d r
e 1 2 e 1 22π2 2π2× κq 1 ; q 2 ; zδμa r;

(13)

where
κq 1 ; q 2 ; z ∞ XlXl 0 m−l

−1m k l m q1 ; zkl m q2 ; L − z: (14)(6)
SOURCESAŝ; ŝ 0 0 Gr; ŝ 0 0 ; r 0 ; ŝ 0 d2 s 0 0  δr − r 0 δŝ − ŝ 0 

(12)

(7)

and obeys the boundary condition (2). Equation (5) is the analog of the Lippmann–Schwinger equation for the RTE. It canbe linearized by the Rytov approximation [1]:I r; ŝ− ln
I 0 r; ŝZ1
d3 r 0 d2 s 0 Gr; ŝ; r 0 ; ŝ 0 δμa r 0 I 0 r 0 ; ŝ 0 :
(8)
I 0 r; ŝ

DETECTORS

y

x

z
0

Fig. 1. Illustrating the slab geometry.

LVol. 33, No. 4 / April 2016 / Journal of the Optical Society of America A

Research Article

Z

3. INVERSE PROBLEM
The inverse problem consists of recovering the absorption δμafrom measurements of the data function ϕ. To proceed, we
assume that the sources and detectors are placed on square lattices with lattice spacings hs and hd , respectively. The positionsof the sources and detectors are given byρs  hs nsx x̂  nsy ŷ;ρd  hd nd x x̂  nd y ŷ;(15)

where nsx , nsy , nd x , and nd y are integers. We introduce theFourier transform of the data function:Xe iqs ·ρs qd ·ρd  ϕρs ; ρd :(16)
ϕ̃q s ; q d  ρs ;ρdδμa ρ; z 

d2 q −iq·ρ X Rσ ν q X
hf qjpi2e
σ 2ν q p;p 0 ν
FBZ 2πν

× K  q; p 0 ; zhp 0 jf ν qiΦq; p:(25)Here, f ν and σ ν are the singular functions and singular valuesof the matrix M ,X
hpjM qjp 0 ihp 0 jf ν qi  σ 2ν qhpjf ν qi;(26)
p0

where M is defined by
Z L
K q; p; zK  q; p 0 ; zdz:
hpjM qjp 0 i 

To compute the Fourier transform, we make use of Eq. (13)and the identity 2 XX
2π
e iq·ρ δq  v;
(18)
h
ρ
v
where v denotes a reciprocal lattice vectorq ∈ − πh ; πh × − πh ; πh. We thus obtainZ
1 X
dzκq s  vs ; −q d  vd ; zϕ̃q s ; q d   2 2
hs hd v ;vanddfa q  q  vs − vd ; z;
× δμ
s
d

(19)

The regularizer R is introduced to limit the effect ofsmall singular values. Common choices of R include the
Tikhonov regularizer Rσ  σ 2 ∕σ 2  ϵ or the step functionRσ  θσ − ϵ, where ϵ > 0. In our computations, both
approaches give similar results.
4. NUMERICAL SIMULATIONSIn this section we report the results of numerical reconstructions from simulated forward data. The medium to be reconstructed consists of a set of five point absorbers, as shown in
Fig. 2. The absorbers are assumed to have the same scattering
coefficients and phase functions as the medium in which they
are embedded, but have different absorption coefficients. The
absorption coefficient δμa is given byδμa r  l2

5Xδr − ri ;

(28)

i1

where
fa q; z 
δμ

Z
d ρe
2

iq·ρ

δμa ρ; z:(20)

Next, we perform the change of variables,qqq d  − p;(21)
q s   p;22
where q s and q d are two-dimensional vectors. We also
assume that δμa is transversely band-limited to the FBZs of
the source-detector lattices (corresponding to the band limit
minf2π∕hs ; 2π∕hd g) and therefore put vs  0 and vd  0
in Eq. (19), which thus becomesZ L
fa q; zdz;
Φq; p K q; p; zδμ(22)

where the positions of the absorbers are x 1 ; y 1 ; z 1  x 2 ; y2 ; z 2   −2l ; 4l ; z 0 − l ,
−4l ; 2l ; z 0 − l ,
x 3 ; y 3 ; z 3   0; 0; z 0 , x 4 ; y 4 ; z 4   2l ; −4l ; z 0  l ,and x 5 ; y 5 ; z 5   4l ; −2l ; z 0  l . The optical properties of the background medium are taken to beμ̄a  0.05 cm−1 ;μs  100 cm−1 ;g  0.9: (29)Note that l  1.00 mm in this case. The above values are
typical for biological tissues in the near-IR spectral range. The
data function Φ is calculated according to Eq. (22) for a sample
with prescribed δμa . Gaussian noise with zero mean and a standard deviation of 1% of the average signal is added to the data
function Φ for each numerical experiment. Reconstructions of

0

where

(27)

0

Here, the vectors q s and q d lie in the first Brillouin zone (FBZ)
of the source and detector lattices:
π
π
π
π
y
y
;
− ≤ qx
:(17)− ≤ q xs ; qs ≤d ; qd ≤
hs
hs
hd
hd

s

553



q
q
Φq; p  h2s h2d ϕ̃  p; − p ;
2
2

(23)


q
q
K q; p; z  κ p  ; p − ; z :
2
2

(24)



For fixed q, Eq. (22) defines a system of one-dimensional
fa . Following the
integral equations for the Fourier transform δμgeneral approach of [23], we construct the pseudoinverse solution of Eq. (22) and perform an inverse Fourier transform toobtain the inversion formula:Fig. 2. Illustrating the set of five point absorbers.

554

Vol. 33, No. 4 / April 2016 / Journal of the Optical Society of America A

δμa are calculated from Eq. (25) using the step function
regularizer. The sum over l in the Green’s function (12) is truncated at l  l max  9.
The source and detector positions are given by Eq. (15) with
−N s ≤ nsx ; nsy ≤ N s and −N d ≤ nd x ; nd y ≤ N d . We also set
N d  120 and N s  120hd ∕hs . The wave vectors q and p
are discretized asq  hd 2N2πFBZ 1 jx x̂  jy ŷ;−N q ≤ jx ; jy ≤ N q ;p  hd 2N2πFBZ 1 kx x̂  ky ŷ;0 ≤ kx ; k y ≤ N p − 1;(30)

where N FBZ is the number of points in the FBZ. We takeN q and N p to be much smaller than N FBZ to control numericalstability at high frequencies. We have found that N FBZ  120,
N q  28, and N p  7 are suitable choices. Finally, we notethat numerical integration over z in Eq. (27) is carried out
using Simpson’s rule, with discretionzation z j  jL − 2l ∕
N z  l j  0; 1; …; N z . For L  10l and 6l we takeN z to be 160 and 80, respectively.

Research Article

the point spread function of the central absorber. Thus the resolution is approximately 2.0l in (a) and 1.5l in (b).
B. Computational Complexity

The reconstructed images are obtained in 90 s using a four-cpu
parallel machine of 1.3 GHz Itanium-2 processors. The mostcomputationally expensive part of the reconstruction is the calculation of the kernel K q; p; z, which is obtained from kl m .
To compute kl m , the Wigner function d lmM and f M n must be
calculated. For each q0 < q ≤ qmax , the construction of d lmM
requires Ol max  1l max  2l max  3 floating point operations and the calculation of f M n requires Ol max  16 
floating point operations. Therefore, when q is approximatedby nq discrete points, the computation of kl m requiresON z nq l 8max  floating point operations. Thus the computationof K requires ON 2q N 2p N z l 2max  floating point operations.
5. EXPERIMENTAL RESULTS

A. Results

A. Experimental Setup

Reconstructions of δμa are shown in Figs. 3 and 4. Results areshown for two different slab thicknesses L  6l , 10l . The
source and detector lattice spacings are hs  0.2l andhd  0.1l , respectively. In Fig. 3, we take L  10l andpanel
z 0  5l in panel (a) and L  6l and z 0  3l inpﬃﬃﬃﬃﬃ(b). The regularization
parameter
is
set
to
be
ε∕
l

pﬃﬃﬃﬃﬃ
−6−5
10 for (a) and ε∕ l  10 for (b). We see that theabsorber is clearly reconstructed in both panels. We will define
the resolution as the full width at half-maximum (FWHM) of

We have tested the reconstruction algorithm with experimental
data using a previously described noncontact OT system [4].
The source is a continuous-wave stabilized diode laser operating
at a wavelength of 785 nm with an output power of 6 mW.
The sample chamber is a rectangular box, constructed of clear
acrylic, of depth 1 cm with square faces of area 15 cm × 15 cm.
The chamber is filled with a scattering medium in which theobjects to be imaged are placed. The medium consists of a 1%
solution of Intralipid in water. The beam is scanned by a pair ofmirrors on one face of the sample and the transmitted beam isimaged onto a CCD.

Fig. 3. Reconstructed δμa x; y  0; z  z 0 ∕ maxδμa x; 0; z 0 
are shown as a function of x∕l . The following cases are considered:L  10l and z 0  5l (top) and L  6l and z 0  3l (bottom).

Fig. 4. Tomographic reconstruction of δμa for five point absorbers
in a slab with L  6l and for z 0  L∕2. The function shown by thecolor scale is δμa r∕ maxr δμa r. The field of view in each panel is
16l × 16l .

Vol. 33, No. 4 / April 2016 / Journal of the Optical Society of America A

Research Article

555

A tomographic dataset is acquired by raster scanning thebeam over a 29 × 29 square lattice with a lattice spacingof 3.21 mm. This yields 841 source positions within a
9 cm × 9 cm area centered on the optical axis. For each source,
a 397 × 397 pixel region of interest is read out from the CCD.
This results in 157,609 detectors arranged in a square latticewith an effective lattice spacing equivalent to 0.36 mm. Thus
a dataset of 1.3 × 108 source-detector pairs is acquired.
B. Two Rods

In this section, we show the reconstruction of a pair of black
metal rods. The rods have a diameter of 3 mm and are suspended in the midplane of the sample chamber. The opticalproperties of the background medium areμ̄a  0.02 cm−1 ;μs  20 cm−1 ;g  0.65;l  20hd ;(31)k  0.02h−1
d ;pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
where k  cμa ∕D and hd  0.36 mm is the lattice spacing
on the detector plane. These parameters yield a transport mean
free path of l  1.4 mm.
Initially, we set the separation of the rods to be 3 mm andtake l max  9. In Fig. 5, tomographic images are displayed witha slice separation of 1 mm. It can be seen that in the central
slice, which is equidistant from the source and the detector
planes, that the rods are well resolved. The shallower anddeeper slices show that the rods remain well resolved but witha smaller diameter, as expected. Figure 6 shows a plot of thereconstruction along the line passing through the centers ofboth rods in the central slice. The distance between the peaks
is 5.8 mm, which slightly underestimates the center-to-center
separation of the rods. The FWHM of the peaks is 2.3 mm,
which underestimates the diameter of the rods.

Fig. 6. Reconstructed δμa ∕μ̄a of the metal rods with (top) 3 mm
spacing and (bottom) 2 mm spacing along the x-axis in the central
plane at y  0 and z  5 mm.

In a second experiment, we set the separation between therods to be 2 mm. In this case the reconstructed image is poorly
resolved. The lower panel in Fig. 6 shows a two-dimensional
plot in the central slice. It is important to note that the reconstructed contrast is not expected to be quantitative owing to thestrong absorption in the interior of the metal rods. However, as
previously noted, the shape of the rods is recovered well.
C. Lemon and Lotus Root Slices

In a final set of experiments, a lemon or lotus root slice of thickness approximately 4 mm is placed in the center of the sample
chamber. The background medium consists of 0.7% Intralipid
plus 0.07% India ink, with optical properties

Fig. 5. Reconstructions of two rods with 3 mm spacing. The fieldof view is 6 cm × 6 cm. The number above each panel shows the
distance from the source plane.

μ̄a  0.4 cm−1 ;μs  20 cm−1 ;k  0.13h−1d ;l  20hd :g  0.65;(32)

Fig. 7. RTE reconstructions of a lemon slice. The field of view is
9 cm × 9 cm. The number above each panel shows the distance fromthe source plane.

556

Vol. 33, No. 4 / April 2016 / Journal of the Optical Society of America A

Research Article

with reflecting boundaries. We note that the MRRF has also
been extended to the case of reflecting boundaries in the halfspace geometry [28].
We begin by noting that the Green’s function for the RTE in
the slab geometry satisfies the equationŝ · ∇  μt Gr; ŝ; ρ0 ; 0; ẑZ μs Aŝ; ŝ 0 Gr; ŝ 0 ; ρ0 ; 0; ẑd2 s 0  δρ − ρ0 δzδŝ − ẑ;Fig. 8. RTE reconstructions of a lotus root slice. The field of view is
9 cm × 9 cm. The number above each panel shows the distance fromthe source plane.

(A1)
together with the boundary conditionsGr; ŝr ; ρ0 ; 0; ẑ  Rjŝi · ẑjGr; ŝi ; ρ0 ; 0; ẑ;

Reconstructions of the lemon and lotus root are shown in
Figs. 7 and 8, respectively. A number of structural features ofthe lemon, including the endocarp, mesocarp, central column
and septa, are readily visualized. Likewise, root tubes of varying
diameter are observable.
6. DISCUSSIONIn summary, we have investigated the inverse problem of
optical tomography in the radiative transport regime. We havedeveloped and tested a fast image reconstruction algorithm thatis applicable to the large datasets of noncontact optical tomography. Our results are illustrated by numerical simulations andreconstructions from experimental data in model systems.
We conclude with several remarks. First, the spatial resolution of the reconstructed images in numerical simulations isclearly superior to that obtained in experiments. There are anumber of contributing factors that explain this finding.
These include the finite size of the source and detector grids,
the presence of shot noise in the detected light, and noise in the
CCD array. We have also not accounted for systematic errors
associated with nonidealities in the optical system, such asreflections from the lens surface. Second, the reconstructions
we have performed are carried out within a linearization of
the RTE. It may be expected that a nonlinear reconstruction
may overcome some of the above limitations, particularly withrespect to quantitative recovery of the absorption. One possible
approach to this problem is inversion of the Born series, which
has been carried out for the case of the DA [20–22]. Inversionof the linearized forward problem is the first step in this procedure, which can be carried out by the method developed
herein. Finally, detailed comparisons of our results with reconstructions carried out within the diffusion approximation willbe the subject of future research. It will be important to understand the role of systematic errors as well as experimental noisein controlling the resolution of reconstructed images. Some discussion of these points has been provided in [4,5].

n̂ · ŝ < 0;z  0; L:(A2)

The reflection coefficient R is given by [29–34]
( h


i
x−nx 0 2
x 0 −nx 2
1
x ≥ xc ;
xnx
nx
2
x
0
0
(A3)
Rx 
1x < xc ;pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
where x 0  1 − n2 1 − x 2 , x c  n2 − 1∕n, and n is theindex of refraction in the slab walls (see Fig. 9). We assume
that the phase function A is given byAŝ; ŝ 0  

l max XlX
l 0 m−l

Al Y l m ŝY l m ŝ 0 ;

where A0  1, A1  g ∈ 0; 1, and Al ∈ 0; 1 (l  2; 3; …).
We adopt the Henyey–Greenstein model [35], and set
Al  g l ;l  0; 1; 2; …:(A5)

In the MRRF [25–28], the specific intensity is given as aand I −superposition of eigenmodes I μμ , whereiq·ρ−Q μ qzI μ r; ŝ; q  el max XlX
−1m −imφq̂
hl jϕμ i
pﬃﬃﬃﬃ eσll 0 m−l

× d lmM iτqλμ Y l m ŝ;

iq·ρQ μ qzI −μ r; ŝ; q  e(A6)l max Xl
X
−1l −imφq̂
hl jϕμ i
pﬃﬃﬃﬃ eσll 0 m−l

× d lm;−M iτqλμ Y l m ŝ;

APPENDIX A: GREEN’S FUNCTION
FOR THE RTE
In this appendix, we outline the calculation of the Green’s function for the RTE using the method of rotated reference frames(MRRF). The MRRF was described in [25–27]. Here we summarize the necessary results for the case of the slab geometry(A4)

Fig. 9. Illustrating the reflecting boundary condition.

(A7)

Vol. 33, No. 4 / April 2016 / Journal of the Optical Society of America A

Research Article
where

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Q μ q  q2  1∕λ2μ :(A8)Here, λn M  and jϕn M i are eigenvalues and eigenvectors ofmatrix BM  defined bysﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
l2 − M2δ0
BM l l 0 2
4l − 1σ l −1 σ l l ;l −1
sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
l  12 − M 2
δ0 ;

(A9)
4l  12 − 1σ l σ l 1 l ;l 1where σ l  μ̄a  μs 1 − g l , l , l 0  jM j; jM j  1; …; l max ,and M  0; 1; 2; …. Note that we make use of amulti-index μ  M ; n, which runs over the setfM ; njλn M  > 0g. The d lmM are analytically continued
Wigner’s d -functions, which can be recursively computed [25].
The Green’s function can be expressed asGr; ŝ; ρ0 ; 0; ẑX Z d2 q−iq·ρ0 
f μ q; nI μ r; ŝ; q
2e
2π
μ
−Q μ qL −
I μ r; ŝ; q:
 f −μ q; ne

(A10)

−
The coefficients f μ q; n and f μ q; n are determinedby the boundary conditions and are solutions to the linear
system   
f q; nq M− q
M
v 

;
M− q M−− q
f − q; n
v −

(A11)
where M ≥ 0, m ≥ 0, l  m  1; m  3; …,
rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
X
2l 0  1
0v l m  δm0 Bl l 0 ∞;v −l m  0;
4π
l0
and
M
l m;M n 

(A12)

X 0
M m
c ll m;M n ; M−−M
l m;M n  −1
l m;M n ; (A13)l0

M m M− ;M−l m;M n  −1
l m;M nX00−−Q μ qL
Ml m;M n  e−1l c ll m;M n ;

(A14)

l0

1
0
l 0 m m
Bl l 0 n pﬃﬃﬃﬃﬃﬃ hl 0 jϕμ ic ll m;M n  Bm
l l 0 ∞ − −1
σl 0
0

0

× fd lmM iτqλμ   1 − δM 0 −1M d lm−M iτqλμ g:
(A15)Here, Bml l 0 n is given by1
Bm
l l 0 n 
2
×sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2l  12l 0  1l − m!l 0 − m!
l  m!l 0  m!

Z 1
0

m
dxP m
l xP l 0 xR n x:

(A16)

(A17)

557

m
m
We note that B−ml l 0 n  Bl l 0 n, and B l l 0 ∞ is definedmby Bl l 0 n with R n x ≡ 1. We also note that f −Mn q; n −1M f M n q; n. Finally, the Green’s function in the slab
geometry is obtained as [25]Gr; ŝ; ρ0 ; 0; ẑZl max Xld2 q −iq·ρ−ρ0  Xe
Y l m ŝi m k l m q; z;≃2π2l 0 m−l(A18)where the functions kl m q; z are given by
X hl jϕn M ikl m q; z  −im e −imφq
pﬃﬃﬃﬃ × e −Q M n qz f M n q; nσl
M ≥0;n −1l mM e −Q M n qL−z f −
M n q; n× d lmM iτqλn M  1 − δM 0 −1M d lm;−M iτqλn M ;(A19)

with φq the angle of q.
Funding. National Institutes of Health (NIH) (P41RR002305, R01EB004832); National Science Foundation(NSF) (DMS-0554100, DMS-1108969).
REFERENCES
1. S. R. Arridge and J. C. Schotland, “Optical tomography: forward andinverse problems,” Inverse Probl. 25, 123010 (2009).
2. G. Turner, G. Zacharakis, A. Soubret, J. Ripoll, and V. Ntziachristos,“Complete-angle projection diffuse optical tomography by use of earlyphotons,” Opt. Lett. 30, 409–411 (2005).
3. R. Schulz, J. Ripoll, and V. Ntziachristos, “Noncontact optical tomography of turbid media,” Opt. Lett. 28, 1701–1703 (2003).
4. Z.-M. Wang, G. Panasyuk, V. Markel, and J. C. Schotland,“Experimental demonstration of an analytic method for image
reconstruction in optical tomography with large data sets,” Opt.
Lett. 30, 3338–3340 (2005).
5. S. D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, andJ. C. Schotland, “Imaging complex structures with diffuse light,” Opt.
Express 16, 5048–5060 (2008).
6. S. Konecky, A. Mahzar, D. Cuccia, A. Durkin, J. Schotland, and B.
Tromberg, “Quantitative optical tomography of subsurface heterogeneities using spatially modulated structured light,” Opt. Express
17, 14780–14790 (2009).
7. G. Bal, “Inverse transport theory and applications,” Inverse Probl. 25,053001 (2009).
8. G. S. Abdoulaev, K. Ren, and A. H. Hielscher, “Optical tomography
as a PDE-constrained optimization problem,” Inverse Probl. 21,
1507–1530 (2005).
9. H. K. Kim and A. H. Hielscher, “A PDE-constrained SQP algorithm
for optical tomography based on the frequency-domain equation ofradiative transfer,” Inverse Probl. 25, 015010 (2009).
10. K. Ren, G. Bal, and A. H. Hielscher, “Frequency domain optical
tomography based on the equation of radiative transfer,” SIAM J.
Sci. Comput. 28, 1463–1489 (2006).
11. O. Dorn, “A transport-backtransport method for optical tomography,”
Inverse Probl. 14, 1107–1130 (1998).
12. P. Gonzalez-Rodriguez and A. D. Kim, “Reflectance optical tomography in epithelial tissues,” Inverse Probl. 25, 015001 (2009).
13. A. Joshi, J. C. Rasmussen, E. M. Sevick-Muraca, T. A. Wareing, andJ. McGhee, “Radiative transport-based frequency-domain fluorescence tomography,” Phys. Med. Biol. 53, 2069–2088 (2008).
14. A. D. Kim, C. Hayakawa, and V. Venugopalan, “Estimating tissue
optical properties using the Born approximation of the transport
equation,” Opt. Lett. 31, 1088–1090 (2006).

558

Vol. 33, No. 4 / April 2016 / Journal of the Optical Society of America A

15. T. Tarvainen, M. Vaukhonen, and S. R. Arridge, “Gauss-Newton
reconstruction method for optical tomography using the finite element
solution of the radiative transfer equation,” J. Quant. Spectrosc.
Radiat. Transfer 109, 2767–2778 (2008).
16. J. Schotland, “Continuous-wave diffusion imaging,” J. Opt. Soc. Am. A
14, 275–279 (1997).
17. J. Schotland and V. Markel, “Inverse scattering with diffusing waves,”
J. Opt. Soc. Am. A 18, 2767–2777 (2001).
18. V. Markel and J. Schotland, “Effects of sampling and limited
data in optical tomography,” Appl. Phys. Lett. 81, 1180–1182 (2002).
19. V. Markel, V. Mital, and J. Schotland, “Inverse problem in optical diffusion tomography. III inversion formulas and singular value decomposition,” J. Opt. Soc. Am. A 20, 890–902 (2003).
20. V. Markel, J. O’sullivan, and J. Schotland, “Inverse problem in optical
diffusion tomography. IV nonlinear inversion formulas,” J. Opt. Soc.
Am. A 20, 903–912 (2003).
21. S. Moskow and J. Schotland, “Convergence and stability of the inverse
scattering series for diffuse waves,” Inverse Probl. 24, 065005 (2008).
22. M. Moskow and J. Schotland, “Numerical studies of the inverse bornseries for diffuse waves,” Inverse Probl. 25, 095007 (2009).
23. V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas,and image reconstruction for optical tomography,” Phys. Rev. E 70,
056616 (2004).
24. J. C. Schotland and V. A. Markel, “Fourier-Laplace structure of theinverse scattering problem for the radiative transport equation,” Inv.
Probl. Imaging 1, 181–188 (2007).
25. M. Machida, G. Y. Panasyuk, J. C. Schotland, and V. A. Markel, “TheGreen’s function for the radiative transport equation in the slab geometry,” J. Phys. A 43, 065402 (2010).

Research Article

26. V. A. Markel, “Modified spherical harmonics method for solving theradiative transport equation,” Waves Random Medium 14, L13–L19(2004).
27. G. Panasyuk, J. C. Schotland, and V. A. Markel, “Radiative transport
equation in rotated reference frames,” J. Phys. A 39, 115–137(2006).
28. A. Liemert and A. Kienle, “Light transport in three-dimensionalsemi-infinite scattering media,” J. Opt. Soc. Am. A 29, 1475–1481(2012).
29. A. D. Kim, “A boundary integral method to compute Green’s functions
for the radiative transport equation,” Waves Random Complex
Medium 15, 17–42 (2005).
30. J. S. Cassell and M. M. R. Williams, “Radiation transport and internalreflection in a sphere,” J. Quant. Spectrosc. Radiat. Transfer 101,
16–28 (2006).
31. J. S. Cassell and M. M. R. Williams, “Radiation transport and internalreflection in a two region, turbid sphere,” J. Quant. Spectrosc. Radiat.
Transfer 104, 400–427 (2007).
32. R. D. M. Garcia, C. E. Siewert, and A. M. Yacout, “On the use of
Fresnel boundary and interface conditions in radiative-transfer calculations for multilayered media,” J. Quant. Spectrosc. Radiat. Transfer
109, 752–769 (2008).
33. R. Aronson, “Boundary conditions for diffusion of light,” J. Opt. Soc.
Am. A 12, 2532–2539 (1995).
34. M. M. R. Williams, “On an integral equation arising in the transport ofradiation through a slab involving internal reflection,” Eur. Phys. J. B
47, 291–293 (2005).
35. L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,”
Astrophys. J. 93, 70–84 (1941).