LUND UNIVERSITY
PO Box 117
221 00 Lund
+46 46-222 00 00
White Monte Carlo for time-resolved photon migration.
Alerstam, Erik; Andersson-Engels, Stefan; Svensson, Tomas
Published in:
Journal of Biomedical Optics
DOI:
10.1117/1.2950319
2008
Link to publication
Citation for published version (APA):
Alerstam, E., Andersson-Engels, S., & Svensson, T. (2008). White Monte Carlo for time-resolved photon
migration.
Journal of Biomedical Optics
,
13
(4), [041304]. https://doi.org/10.1117/1.2950319
Total number of authors:
3
General rights
Unless other specific re-use rights are stated the following general rights apply:
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors
and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the
legal requirements associated with these rights.
• Users may download and print one copy of any publication from the public portal for the purpose of private study
or research.
• You may not further distribute the material or use it for any profit-making activity or commercial gain
• You may freely distribute the URL identifying the publication in the public portal
Read more about Creative commons licenses: https://creativecommons.org/licenses/
Take down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove
access to the work immediately and investigate your claim.
White Monte Carlo for time-resolved photon migration
Erik Alerstam
Stefan Andersson-Engels
Tomas Svensson
Lund University
Department of Physics
Sweden
Abstract. A novel scheme for fully scalable White Monte Carlo
共WMC兲 has been developed and is used as a forward solver in the
evaluation of experimental time-resolved spectroscopy. Previously re-
ported scaling problems are avoided by storing detection events indi-
vidually, turning spatial and temporal binning into post-simulation ac-
tivities. The approach is suitable for modeling of both interstitial and
noninvasive settings 共i.e., infinite and semi-infinite geometries兲. Moti-
vated by an interest in in vivo optical properties of human prostate
tissue, we utilize WMC to explore the low albedo regime of time-
domain photon migration—a regime where the diffusion approxima-
tion of radiative transport theory breaks down, leading to the risk of
overestimating both reduced scattering
共
s
⬘
兲 and absorption 共
a
兲. Ex-
perimental work supports our findings and establishes the advantages
of Monte Carlo–based evaluation.
© 2008 Society of Photo-Optical Instrumenta-
tion Engineers. 关DOI: 10.1117/1.2950319兴
Keywords: Monte Carlo; photon migration; time-resolved spectroscopy 共TRS兲; opti-
cal properties; human prostate
.
Paper 07390SSR received Sep. 18, 2007; revised manuscript received Nov. 13,
2007; accepted for publication Nov. 16, 2007; published online Jul. 9, 2008.
1 Introduction
Many applications within the field of biomedical optics rely
either on the capability of performing, or the availability of,
accurate measurements of optical properties of highly scatter-
ing materials. This is reflected by the massive development of
theory related to light propagation in turbid media 共photon
migration兲, and the numerous techniques available for charac-
terisation of such materials.
Time-resolved spectroscopy 共TRS兲 is one of several tech-
niques available for assessing optical properties of turbid me-
dia. By studying the broadening of short 共picosecond兲 light
pulses, it allows determination of both the absorption coeffi-
cient
共
a
兲 and the reduced scattering coefficient 共
s
⬘
兲 without
the need of absolute measurements of light intensities. At first,
time-resolved photon migration was investigated and evalu-
ated using the diffusion approximation of radiative transport
theory.
1
Although proven useful in numerous cases, the diffu-
sion modeling was found erroneous at low albedos or close to
radiation sources.
2,3
Unfortunately, several tissue types exhibit
optical properties in the range where the use of the diffusion
approximation may be questioned. The human prostate is one
example,
4,5
exhibiting reduced scattering below 10 cm
−1
and
absorption above
0.3 cm
−1
.
In order to extend the range of optical properties and
source-detector separations over which TRS can provide ac-
curate data, methods for Monte Carlo-based modeling were
developed. Introduced to the field of biomedical optics and
photon migration by Wilson and Adam,
6
Monte Carlo 共MC兲
simulation has become the gold standard for modeling of light
propagation in tissue optics.
7
Besides modeling spatially and
temporally resolved light distribution, MC has proven useful
in, for example, fluorescence modeling
8
and Raman
spectroscopy.
9
Traditionally, MC was performed for a particular set of
optical properties at a time. Since the simulation is time-
consuming, it is thus not useful as a forward solver in reverese
problems, for instance, during evaluation of experimental data
共iterative curve-fitting兲. This obstacle lead to the development
of White Monte Carlo 共WMC兲,
10–12
in which a single simula-
tion in combination with proper rescaling ensures coverage of
a wide range of optical properties. The main feature in WMC,
making it feasible for use as a forward model in an iterative
solver, is illustrated in Fig. 1.
During MC simulations of light propagation in homoge-
neously scattering and nonabsorbing media, photon paths are
determined only by the scattering phase-function, the scatter-
ing coefficient, and the sequence of random numbers gener-
ated by the simulation program. For a given phase function,
and considering a particular sequence of random numbers, the
photon path will scale linearly with the scattering coefficient,
s
. The recording of photon paths and corresponding time-of-
flights thus allows post-simulation transformation, from the
scattering coefficient used during simulation, to an arbitrary.
Furthermore, the impact of nonzero absorption can be im-
posed post-simulation simply by giving photons different
weights
w
a
according to the Beer-Lambert law of attenuation.
This weight is stated in Eq. 共1兲, where
c
⬘
is the speed of light
within the media:
w
a
= exp共−
a
c
⬘
t兲. 共1兲
This means that a single Monte Carlo simulation can be used
to extract photon time-of-flight distribution not only for dif-
1083-3668/2008/13共4兲/041304/10/$25.00 © 2008 SPIE
Address all correspondence to: Tomas Svensson, Lund University, Department
of Physics, Lund, Sweden.
Journal of Biomedical Optics 13共4兲, 041304 共July/August 2008兲
Journal of Biomedical Optics July/August 2008
쎲
Vol. 13共4兲041304-1
Downloaded from SPIE Digital Library on 01 Jul 2011 to 130.235.188.41. Terms of Use: http://spiedl.org/terms
ferent source-detector separations, but also for different opti-
cal properties
s
⬘
and
a
.
The preceding theory has been discussed in several
papers.
8,10–15
Graaff et al. suggested a limited scalable Monte
Carlo technique where the optical properties of simulations in
slab geometries could be scaled as long as the total attenua-
tion coefficient was held constant.
13
Two groups simulta-
neously and independently extended the theory of Graaff et al.
Kienle and Patterson suggested a scalable Monte Carlo tech-
nique for infinite and semi-infinite homogeneously scattering
media and used independent 共reference兲 MC simulations for
verification of its performance in the semi-infinite case.
10
Pifferi et al. suggested a similar approach,
11,12
proposing scal-
ability in both
a
and
s
⬘
. In practice, however, they based
their evaluation on interpolation between MC simulations car-
ried out at different
s
, thus utilizing scalability in
a
only.
On the other hand, their work included evaluation of experi-
mental time-resolved transmittance in the slab geometry 共a
geometry not allowing
s
rescaling兲. Swartling et al. showed
the usefulness of the WMC approach in fluorescence emission
spectra modeling.
8
Swartling also raised the important ques-
tion regarding the equivalency of WMC and traditional MC.
14
Xu et al. demonstrated the superior performance of WMC
over different light propagation models as a forward model
for evaluation of frequency domain data 共generated by a tra-
ditional Monte Carlo program
15
兲. In the same paper, Xu et al.
demonstrated scaling of data from an absorbing media using
the weighting relationships from traditional Monte Carlo. Xu
et al. also reported the so-called scaling effect as being the
major inconvenience in WMC-based data evaluation. This in-
convenience originates from the fact that scaling in
s
is ac-
companied by a scaling of temporal and spatial bin size, re-
sulting in a need for data resampling and a limited range in
scalability.
Motivated by an interest in in vivo optical characterization
of human prostate tissue, this work is aimed at providing a
scheme for fully scalable WMC for time-domain photon mi-
gration and demonstrating its value in evaluation of experi-
mental data in the low albedo regime of photon migration.
The approach is useful in both infinite and semi-infinite ge-
ometries, featuring individual storage of the spatial location
and the time-of-flight for all potential detection events. Since
this allows post-scaling binning 共temporal, as well as spatial兲,
it eliminates the need for the data resampling otherwise ac-
companying
s
scaling. It also allows an accurate account for
finite extension of source and detector areas 共e.g., optical fiber
diameters兲.
2 Materials and Methods
2.1 White Monte Carlo
A WMC model was developed to serve as the forward model
for evaluation of fiber-based time-resolved spectroscopy un-
der interstitial as well as under noninvasive conditions 共i.e.,
infinite and semi-infinite geometry兲. The model consists of a
simulation program written in C and a set of
MATLAB scripts
performing post-simulation processing. The main objectives
were to retain full scalability in both
s
⬘
and
a
, while avoid-
ing scaling inconveniences by moving spatial and temporal
binning to post-simulation. This eliminates the need for any
temporal resampling and allows accurate convolutions to
properly account for the finite size 共radius
R
f
兲 and light dis-
tribution of optical fibers. A schematic illustration of the
WMC model is given in Fig. 2. Apart from selecting simula-
tion geometry, the user has to provide the WMC simulation
program with a few input parameters. The refractive index
n,
the anisotropy factor
g, and the numerical aperture NA of the
involved optical fibers are material parameters specific to a
µ
s
/α
αtt
αrr
µ
s
Fig. 1 In WMC, the shapes of the photon paths are determined only
by the phase scattering function and the sequence of random num-
bers. Paths generated by simulations in nonabsorbing media can
hence be linearly scaled to apply for different values of
s
. The illus-
tration is adapted from Ref. 12.
WMC Simulation
Sorting
Database range:
Fibre properties:
POST SIMULATION
FORWARD MODELUSER-SUPPLIED DATA
WMC INPUT DATA
SIMULATION
Scale & select
Calculate
geometric weights
Calculate
absorption weights
Temporal binning
Time-of-flight
histogram
∆t
R
f
R
f
ρ
µ
a
µ
s
ρ
g
and n
NA
Material constants:
µ
max
s
t
max
and
Fig. 2 A flowchart of the WMC scheme used. The WMC simulation
program performs the simulation in either infinite or semi-infinite ho-
mogenously scattering media, using user-supplied WMC input data.
The resulting database is sorted and stored to be used later by the
forward model. The forward model rescales the database data and
generates time-of-flight histograms corresponding to parameters sup-
plied by the user. These parameters include the radius R
f
of the in-
volved optical fibers, the source-detector separation
, the optical
properties
s
⬘
and
a
, and the desired temporal channel width ⌬t.As
this fast forward model is incorporated in the forward solver, the
method can be used to evaluate experimental time-domain data.
Alerstam, Andersson-Engels, and Svensson: White Monte Carlo for time-resolved photon migration
Journal of Biomedical Optics July/August 2008
쎲
Vol. 13共4兲041304-2
Downloaded from SPIE Digital Library on 01 Jul 2011 to 130.235.188.41. Terms of Use: http://spiedl.org/terms
particular simulation. The simulation is run at the specified
scattering level
s
max
but is always terminated when time
reaches
t
max
. This notation is used because
s
max
defines the
maximum scattering coefficient for which the resulting data-
base can provide valid data throughout the time interval
关0,t
max
兴. Generally, scaling to a certain scattering
s
=
s
max
/
␣
results in data valid in the time interval 关0,
␣
t
max
兴.
2.1.1 Simulation program
The WMC simulation program was written in ANSI C, adapt-
ing some parts from the de facto standard MC simulation
program multi-layer Monte Carlo 共MCML
7
兲. As the photon
propagation is similar to traditional MC, we refer to the work
of Wang et al.
7
and Prahl et al.
16
for basic Monte Carlo theory.
Here, we present features differing from their work.
An important change is the adaptation of a state-of-the-art
pseudo-random number generator, the Mersenne twister by
Saito and Matsumoto.
17
The implementation used was the
double-precision SIMD-oriented fast Mersenne twister
共dSIMD兲 version 1.2.1, featuring a
2
132049
−1 period, docu-
mented excellent equidistribution properties, and fast random
number generation. The entire
32-bit output of the time ()
function in C was used to seed the generator.
Photons are launched from the origin of a Cartesian coor-
dinate system 共as in MCML兲. The launch direction is in the
positive
z direction 共downward兲 with the addition of a deflec-
tion angle,
, representing the emission cone of the source
fiber, defined by the NA of the fiber,
max
=arcsin共NA/ n兲. The
angular distribution is assumed flat, i.e.,
cos
=1−
共1
−cos
max
兲, where
is a random number,
⑀
关0,1兲. Although
unnecessary in a cylindrically symmetric problem, the azi-
muthal angle
is also randomized,
=2
, where
⑀
关0,1兲.
Note that taking the angular distribution of the source fiber
into account prevents post-simulation scaling of the refractive
index. Thus, different WMC simulations are required in order
to handle different fiber
NA, as well as different refractive
indices.
The step size,
⌬s, is calculated using the scattering coeffi-
cient instead of the total attenuation coefficient. It is defined
in Eq. 共2兲:
⌬s =
−ln
共
兲
s
max
, 共2兲
where
is a random variable,
⑀
共0,1兴, implying ⌬s
⑀
关0,⬁兲.
As in MCML, the actual photon scattering is simulated using
the Henyey-Greenstein phase function.
The photon detection scheme assumes that the source and
detector optical fibers are parallel, having their tips in the
same plane 共a plane to which the fibers are also assumed to be
orthogonal兲. This corresponds to the settings followed in the
interstitial clinical measurements on human prostate tissue
presented in Ref. 5, where optical fibers are inserted to the
same depth. Regardless of whether simulations are performed
in the infinite or semi-infinite geometry, a potential photon
detection event requires that a photon path crosses the source-
detector plane.
In the case of an infinite medium, photons passing upward
through the source-detector plane, being within the accep-
tance cone of the detector fiber, may be detected. Such a
crossing is referred to as a detection event and is always reg-
istered by storing its radial distance from the source
r, as well
as its total time-of-flight
t. For memory conservation, storing
is done using
32-bit floating point variables. Since the photon
is propagating in an infinite medium, and since the position of
the detector fiber is undefined, the photon is, however, not
terminated at this point. Instead, photon termination occurs
only when the time reaches
t
max
. Note that this implies that a
single photon may generate multiple detection events. These
events are considered independent by the proposed WMC
scheme, inducing a small error in generated time-of-flight his-
tograms. This issue is further discussed in Sec. 2.1.3.
For the case of a semi-infinite geometry, the source-
detector plane equals the medium boundary. Photons escaping
the medium at an angle within the acceptance cone of the
detector fiber generate detection events, i.e.,
r and t are re-
corded. Here, a notable difference from traditional Monte
Carlo in semi-infinite geometries is that partial reflections are
not allowed by WMC 共the weight of the photons are not
monitored兲. Instead, the reflection coefficient is compared to a
random number, and the photons are either reflected or trans-
mitted and terminated. Note that this also means that the
semi-infinite detection geometry does not suffer from the risk
of double-detecting photons. Note also that photons are termi-
nated not only upon boundary crossing, but also when the
time exceeds the predefined maximum time-of-flight
t
max
.
2.1.2 Post-simulation processing
The database of detection events 共r, t pairs兲 from the simula-
tion are sorted with respect to
r. This is done using the fast
兵
O关n log共n兲兴 worst-case time complexity其 in-place sorting al-
gorithm Combsort11.
18
The O共 1兲 memory usage of the
Combsort11 algorithm enables sorting and usage of databases
of sizes close to the physical memory of the computer used.
The sorted database is used by a
MATLAB function gener-
ating photon time-of-flight distributions. This function is em-
ployed as the forward model in an iterative solver and is
called with the following input parameters:
a
and
s
⬘
being
the optical properties corresponding to the resulting time-of-
flight distribution, the temporal channel width
⌬t, source-
detector separation
, and the fiber radius R
f
共experiments
typically involve two identical fibers兲. The first task is to cal-
culate the scaling coefficient
␣
=
s
max
/
s
共
␣
艌1兲. The r, t
pairs within the spatial interval
−2R
f
⬍
␣
r⬍
+2R
f
are ex-
tracted and scaled 共
r
⬘
=
␣
r, t
⬘
=
␣
t兲. As the r, t pairs are sorted
with respect to
r, the selecting process turns into a simple task
of finding the borders of the spatial interval, which is done
quickly using a standard binary search algorithm. The interval
itself is motivated by the convolution of each detection event
with the size of the source fiber combined with the size of the
detector fiber. A photon will hence contribute to the time dis-
persion histogram as soon as its distance to the center of the
detector fiber is less than
2R
f
. Each photon detection event is
assigned a weight,
w
f
共r
⬘
兲, according to its spatial position.
This weight is based on the overlap between two circles of
radius
R
f
, one centered at the detection event and one cen-
tered at the detector location. More precisely, the weight is
reached after integration over all origin-to-detector angles 共ro-
tating the detector circle around the
z axis兲. This is similar to
work by Wang et al.
19
and Prahl.
20
As the fiber radius and the
Alerstam, Andersson-Engels, and Svensson: White Monte Carlo for time-resolved photon migration
Journal of Biomedical Optics July/August 2008
쎲
Vol. 13共4兲041304-3
Downloaded from SPIE Digital Library on 01 Jul 2011 to 130.235.188.41. Terms of Use: http://spiedl.org/terms
source-detector separation is known when evaluating data, the
weight
w
f
共r
⬘
兲 can be precalculated, and the computationally
costly integrals do not have to be evaluated for each photon. A
uniform near-field irradiance distribution of the fibers was as-
sumed in this work, although any distribution could be uti-
lized with minimal modifications.
The photon detection events are also assigned a weight due
to absorption,
w
a
共t
⬘
兲= exp共−
a
c
⬘
t
⬘
兲, where c
⬘
is the local
speed of light. The total weight for each detection event is
thus
w
tot
=w
f
·w
a
.
The final step is a simple matter of creating a weighted
histogram of the 共
t
⬘
, w
tot
兲 data, where the temporal channel
width,
⌬t, can be specified by the user.
2.1.3 Entangled detection events
As mentioned earlier, a single photon may generate multiple
detection events 共a primary event, a secondary event, and so
on兲. All such events should contribute to the final time-of-
flight histogram 共given that they are located within the detec-
tion range
−2R
f
⬍r
⬘
⬍
+2R
f
兲. However, if the spatial dis-
tance between two such events is less than
2R
f
, their weights
cannot be considered independent, and there is a risk for er-
roneous additions of weight to the final histogram. The worst
case is if the two events coincide spatially 共zero event dis-
tance兲, a case in which the entire weight contributed by the
secondary event is erroneous. In contrast, when the spatial
event separation is
2R
f
or greater, none of its contribution is
erroneous. For event separations between zero and
2R
f
, the
matter is nontrivial. In general, the erroneously added weight
decreases with event separation. An estimate of the total
amount of erroneous weight is reached by studying the mul-
tiple detection events occurring during the WMC simulations.
For the case of
s
⬘
=10 cm
−1
, g= 0.7, NA=0.29, and R
f
=0.3 mm, we found that 0.5% of the primary events are ac-
companied by a secondary event less than
2R
f
away 共and only
0.1% within
R
f
兲. Assuming that all corresponding weight is
erroneous 共clearly an exaggeration of the problem兲,wecan
note that the final photon time-of-flight histogram contains
less than 0.5% erroneous weight. Although this suggests that
the problem of entangled detection events is insignificant, it
may deserve further attention.
2.1.4 Verification
In an effort to validate the WMC simulation program code
and the implementation of the scaling relationships, the open
source program MCML was carefully modified to monitor the
time-of-flight histograms of photons leaving a homogenously
scattering, semi-infinite medium. The photon weights leaving
the medium within the
1-mm-wide spatial bins, centered
around 10, 15, and
20 mm, were recorded with 10-ps tempo-
ral resolution.
10
8
photons were simulated at g = 0.7 and n
=1.33
for all combinations of
s
⬘
=5, 10 cm
−1
and
a
=0.1,
0.5 cm
−1
.
2.2 Time-Resolved Instrumentation
Time-resolved photon migration experiments were conducted
using a compact
共50⫻ 50⫻ 30 cm
3
兲 and portable time-
domain photon migration instrument primarily intended for
spectroscopy of biological tissues in clinical environments.
Detailed information on the instrumentation can be found in
previous publications.
5,21,22
Briefly, the system is based on
pulsed diode laser technology and time correlated single pho-
ton counting 共TCSPC兲. Four pulsed diode lasers 共at 660, 786,
830, and
916 nm, respectively兲 generate pulses having a
FWHM of about
70 ps, the average power being 1to2mW.
Light is injected into the sample and collected using
600-
m
GRIN optical fibers 共NA=0.29兲. A fast MCP-PMT together
with a TCSPC computer card is used to obtain photon time-
of-flight histograms. Broadening in fibers and detector yields
an instrument response function 共IRF兲 of about
100-ps
FWHM. The IRF is measured by putting the fiber ends face to
face with a thin paper coated on both sides with black toner,
similar to the IRF measurements proposed by Schmidt et al.
23
A schematic illustration of the setup is given in Fig. 3.
2.3 Experiment
In order to compare diffusion-based and WMC-based evalua-
tion of experimental data, measurements were performed on
tissue phantoms based on Intralipid 共Fresenius Kabi
200 mg/ ml兲 and ink 共1:100 stock solution of Pelikan Fount
India ink兲.
24
Two measurement series were performed: one added ab-
sorber series using ink,
25
and one added scatterer series using
Intralipid.
26,27
Neglecting the minor volume change occurring
during these measurement series, we can qualitatively expect
a constant scattering and linearly increasing absorption in the
added absorber series. In addition, when extrapolating of de-
rived
a
toward the zero ink level, the absorption should be
close to the absorption exhibited by pure water.
Similarly, the added scatterer series should yield a constant
absorption, linearly increasing scattering. Extrapolation of the
derived scattering coefficient toward the zero Intralipid level
should yield a zero scattering 共at least if ink scattering is
negligible兲.
The added absorber series was performed on a phantom
based on
577 ml water and 22 ml Intralipid. Measurements
were performed at
2to8mltotal volume of added ink solu-
tion 共
1-ml increments兲. The added scatterer series was per-
formed on a phantom based on
577 ml water and 4mlink
solution. Measurements were performed at
10 to 24 ml added
volume of Intralipid 共in
2-ml increments兲. A cylindrical plas-
tic container,
Ø= 110-mm and height=70 mm, was used to
hold the phantoms during mixing and measurements. The
container was placed on a magnetic stirrer, utilized for mixing
of the phantoms after each addition of ink or Intralipid. Fibers
were placed in parallel at
30-mm depth in the center of the
phantom, separated
14.7 mm, center to center. Data acquisi-
tion was performed for
30 s for each measurement, and the
IRF was measured before, after, and occasionally between the
measurements 共monitoring potential temporal drifts in the
system兲.
Laser Driver
660 nm
786 nm
830 nm
916 nm
Sample
TCSPC
MCP-PMT
SYNC
Cooling
Fig. 3 A schematic of the instrumentation in interstitial mode.
Alerstam, Andersson-Engels, and Svensson: White Monte Carlo for time-resolved photon migration
Journal of Biomedical Optics July/August 2008
쎲
Vol. 13共4兲041304-4
Downloaded from SPIE Digital Library on 01 Jul 2011 to 130.235.188.41. Terms of Use: http://spiedl.org/terms
