Class 3 - Monte Carlo 1

Report
Modeling of polarized light
transfer in scattering media
Monte Carlo
II Escuela de Optica Biomedica, Puebla, 2011
Introduction
The Monte Carlo method simulates the “random
walk” of photons in a medium that contains
absorption and scattering.
II Escuela de Optica Biomedica, Puebla, 2011
Introduction
• The method is based on a set of rules that
governs the movement of a photon in tissue.
• The two key decisions are:
•
•
The mean free path for a scattering or
absorption event.
The scattering angle.
II Escuela de Optica Biomedica, Puebla, 2011
Introduction
• Monte Carlo programs estimates ensamble average
quantities
• The Monte Carlo method is essentially a technique
for sampling a probability density function based on
a computer generated random number.
• Need large number of photons
– Error )snotohp fo rebmuN(√/1
II Escuela de Optica Biomedica, Puebla, 2011
Monte Carlo flow chart
II Escuela de Optica Biomedica, Puebla, 2011
Polarized light Monte Carlo
• Main issues
– Define a polarization reference plane
– Propagate a Stokes vector & reference plane
– Handle scattering – No Heyney Greenstain phase
function
– Handle scattering – Rejection method/Mie
scatterers
– Handle interface and exit
II Escuela de Optica Biomedica, Puebla, 2011
J.C. Ramella-Roman, S.A. Prahl , S.L. Jacques, 'Three Monte Carlo programs of polarized light
transport into scattering media: part I,' Optics Express, 13, pp. 4420-4438, 2005.
J.C. Ramella-Roman, S.A. Prahl, S.L. Jacques, 'Three Monte Carlo programs of polarized light
transport into scattering media: part II,' Optics Express, 13, pp.10392-10405, 2005.
II Escuela de Optica Biomedica, Puebla, 2011
Three methods
• Meridian Plane Monte Carlo – Today
• Euler and Quaternion Monte Carlo –
Tomorrow
• The difference among three programs is the
tracking of the reference plane
II Escuela de Optica Biomedica, Puebla, 2011
Meridian plane Monte Carlo
II Escuela de Optica Biomedica, Puebla, 2011
Meridian plane
• Chandrasekhar was the first to envision the
scattering of a photon from one location to
the next in the form of meridian and
scattering planes.
• The photon directions, before and after
scattering, are represented as points on the
unit sphere and each photon direction is
uniquely described by two angles q and f.
II Escuela de Optica Biomedica, Puebla, 2011
Meridian plane
• q is the angle between the initial
photon direction and the Z-axis.
– Scattering angle
• f is the angle between the meridian
plane and the X-Z plane.
– Azimuth angle
• The photon direction is specified by
a unit vector I1
• Elements of I1 are the direction
cosines [ux, uy, uz]
II Escuela de Optica Biomedica, Puebla, 2011
Meridian plane
• The unit vector I1 and the Z-axis
determine a plane COA.
• This plane is the meridian plane
II Escuela de Optica Biomedica, Puebla, 2011
Launch
• A photon is launched in a specific
meridian/reference frame and status of
polarization.
• Slab geometry, perpendicular illumination,
pencil beam
• Direction cosines
– [ux, uy, uz] = [0, 0, 1].
II Escuela de Optica Biomedica, Puebla, 2011
Launch-STEP A
• Initially the meridian plane f =0,
• The reference frame is initially equal to the x-z
plane.
• The initial Stokes vector is relative to this
meridian plane.
II Escuela de Optica Biomedica, Puebla, 2011
Launch-STEP B
• The status of polarization as an appropriate
Stokes vector, S = [I Q U V].
• Example S=[1 1 0 0]
– launched photon will be linearly polarized parallel
to the x-z plane.
II Escuela de Optica Biomedica, Puebla, 2011
Move
• The photon is moved as in a standard Monte
Carlo program.
• The photon is moved a propagation distance
Ds that is calculated based on pseudo random
number  generated in the interval (0,1].
Ds  

ln( )
t
µt = µa +µs,
µa is the absorption coefficient of the media.
µs is the scattering coefficient of the media.
II Escuela de Optica Biomedica, Puebla, 2011
Move cnt.
• The trajectory of the photon is specified by
the unit vector I characterized by the direction
cosines [ux, uy, uz].
ux  I x
uy  I y
uz  I z
• where x y and z are the unit vectors in the
laboratory frame XYZ.
II Escuela de Optica Biomedica, Puebla, 2011
Move cnt.
• The photon position is updated to a new
position [x',y',z'] with the following equations:
x x  ux Ds  I x Ds
y y  u y Ds  I y Ds
z z  uz Ds  I z Ds
II Escuela de Optica Biomedica, Puebla, 2011
Scatter
• Three matrix multiplications are necessary to
handle the scattering of a photon and track its
polarization
• 1. Stokes vector rotation in the scattering
plane
• 2. Scattering -updating the Stokes vector
• 3. Return to a new meridian plane
II Escuela de Optica Biomedica, Puebla, 2011
Rotation to a new reference plane
• The E field is originally
defined respect to a
meridian plane COA.
• The field can be
decomposed into its parallel
and perpendicular
components E|| and E.
• The choice of f and q is
done with the rejection
method
II Escuela de Optica Biomedica, Puebla, 2011
Rotation to a new reference plane
• First the Stokes vector is
rotated so that its
reference plane is ABO
• This rotation is necessary
because the scattering
matrix, that defines the
elastic interaction of a
photon with a sphere, is
specified with respect to
the frame of reference of
the scattering plane
II Escuela de Optica Biomedica, Puebla, 2011
1. Stokes vector rotation in the
scattering plane
• The Stokes vector S1
defined relative to the
scattering plane (ABO) is
found by multiplying by a
rotational matrix R(i1)
• This action corresponds
to a counterclockwise
rotation of an angle i1 
about the direction of
propagation.
1
0
0

0 cos( 2i1 ) sin(2i1 )
R i1  
0  sin(2i1 ) cos( 2i1 )

0
0
0

II Escuela de Optica Biomedica, Puebla, 2011
0

0
0

1
Scattering
• The Stokes vector is
multiplied by the
scattering matrix that
accounts for scattering
of the photon at an
angle a
II Escuela de Optica Biomedica, Puebla, 2011
2. Scattering -updating the Stokes
vector
• The Stokes vector is
multiplied by a
scattering matrix, M(a).
s ( a )
11

• a is the scattering angle
s (a )
M a   12
 0
between the direction

 0
of the photon before
and after scattering.
s12 ( a )
s11( a )

II Escuela de Optica Biomedica, Puebla, 2011
0
0


0
0 
s33( a ) s34( a ) 

s34( a ) s33( a ) 
0
0
Updating the Stokes vector
• The parameters s11 , s12
, s33 , s34 are calculated
with Rayleigh or Mie
theory.
• These terms are
expressed as
2
2
1
s11 ( a )  [ S2 ( a )  S1( a ) ]
2
2
2
1
s12 ( a )  [ S2 ( a )  S1( a ) ]
2
1
s33 ( a )  [ S2* ( a )S1( a )  S2 ( a )S1* ( a )]
2
i
s34 ( a )  [ S2* ( a )S1( a )  S2 ( a )S1* ( a )]
2
*http://omlc.ogi.edu/calc/mie_calc.html

• S1, S2* – Mie scattering
C. Bohren and D. R. Huffman,
Absorption and scattering of light by small
particles, (Wiley Science Paperback Series,1998).
II Escuela de Optica Biomedica, Puebla, 2011
Rejection method
• The phase function P(α,β) for incident light
with a Stokes vector So= [Io, Qo, Uo, Vo] is
P(a , b ) = s11 (a ) Io + s12 (a ) [Qo cos (2b ) +Uo sin (2b )]
• α angle of scattering
• b angle of rotation into the scattering plane
II Escuela de Optica Biomedica, Puebla, 2011
Rejection method
• For unpolarized light [1 0 0 0]
P(a , b ) = s11 (a ) Io
• The rejection method is used to generate random

variables with a particular distribution.
II Escuela de Optica Biomedica, Puebla, 2011
Parallel incidence- parallel detection
F
B
D=0.01µm
D=1.0µm
II Escuela de Optica Biomedica, Puebla, 2011
D=2.0µm
Parallel incidence- perpendicular
detection
D=0.01µm
D=1.0µm
II Escuela de Optica Biomedica, Puebla, 2011
D=2.0µm
Rejection method-step 1
• Generate Prand - uniform random number
between 0 and 1.
• arand is generated uniformly between 0 and π.
• The angle arand is accepted as the new scattering
angle if Prand ≤ P(arand).
• If not, a new Prand and arand are generated
II Escuela de Optica Biomedica, Puebla, 2011
Rejection method-step 2
P(a , b ) = s11 (a ) Io + s12 (a ) [Qo cos (2b ) +Uo sin (2b )]
• Generate Prand - uniform random number
between 0 and 1.
• brand is generated uniformly between 0 and 2π.
• The angle brand is accepted as the new scattering
angle if Prand ≤ P(arand, brand).
• If not, a new Prand and arand , brand are generated
II Escuela de Optica Biomedica, Puebla, 2011
Scattering from spherical particles
Single scattering
II Escuela de Optica Biomedica, Puebla, 2011
Update of direction cosines
If |uz| ≈ 1
^
u x  sin(q ) cos( f )
• q and f are obtained at
in the scattering step of
the Monte Carlo
program
^
u y  sin(q ) sin( f )
^
u z  cos(q )
uz
uz
II Escuela de Optica Biomedica, Puebla, 2011
Update of direction cosines
If |uz| ≠ 1
^
ux 
1
1 u2
sin(q )( ux u y cos(f )  u y sin(f ))  ux cos(q )
z
^
uy 
1
1 u2
sin(q )( ux u z cos(f )  ux sin(f ))  u y cos(q )
z
^
u z  1 u2 sin(q ) cos(f )( u y u z cos(f )  ux sin(f ))  u z cos(q )
z
II Escuela de Optica Biomedica, Puebla, 2011
3. Return to a new meridian plane
• The Stokes vector is
rotated so that it is
referenced to the new
meridian plane COB
II Escuela de Optica Biomedica, Puebla, 2011

3. Return to a new meridian plane
• The Stokes vector is
multiplied by the
rotational matrix R(-i2)
so that it is referenced
to the meridian plane
COB.
^
cosi2 
u z  u z cosa
^
 ( 1  cos2 a )( 1  u 2z )
u are the direction cosines before scattering and û are
the direction cosines
after scattering.
II Escuela de Optica Biomedica, Puebla, 2011
3. Return to a new meridian plane
• The Stokes vector is
multiplied by the
rotational matrix R(-i2)
so that it is referenced
to the meridian plane
COB.
1
0
0

0 cos(2i2 ) sin(2i2 )
R i2  
0  sin(2i2 ) cos(2i2 )

0
0
0

0

0
0

1
II Escuela de Optica Biomedica, Puebla, 2011
Summary
• In summary the Stokes vector Snew after a
scattering event is obtained from the Stokes
vector before the scattering event S
Snew  R( i2 )M( a ) R( i1 )S
II Escuela de Optica Biomedica, Puebla, 2011
Choice of f and a
• A fundamental problem in every Monte Carlo
program with polarization information is the
choice of the angles a (angle of scattering)
and f (angle of rotation into the scattering
plane).
• These angles are selected based on the phase
function of the considered scatterers and a
rejection method (**)
** W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in C, the
art of Scientific Computing,IICambridge
University
Press;
(1992).
Escuela de Optica
Biomedica,
Puebla,
2011
Photon life
• The life of a photon ends when the photon
passes through a boundary or when its weight
W value falls below a threshold.
• Roulette is used to terminate the photon
packet when W  Wth.
– Gives the photon packet one chance of surviving
– If the photon packet does not survive the roulette,
the photon weight is reduced to zero and the
photon is terminated.
II Escuela de Optica Biomedica, Puebla, 2011
Boundaries
• One last rotation of the Stokes vector is
necessary to put the photon polarization in
the reference frame of the detector.
• For simplicity we are neglecting Fresnel
reflectance
II Escuela de Optica Biomedica, Puebla, 2011
Boundaries
• For a photon backscattered from a slab, the last
rotation to the meridian plane of the detector is
of an angle y
1
y  tan (
uy
ux
)
• The Stokes vector of the reflected photon is
multiplied one final time by R(y). For a
transmitted photon the angle y
1
y   tan (
uy
ux
)
II Escuela de Optica Biomedica, Puebla, 2011
Results
• Average Stokes vector values
• HI, HQ, HU, HV, VI, VQ )…( IMAGES
– All you need to build an image of the Mueller
matrix
II Escuela de Optica Biomedica, Puebla, 2011
Monte Carlo testing
II Escuela de Optica Biomedica, Puebla, 2011
Comparison with Evans’ addingdoubling*
Diameter
I Evans
I This code
Q Evans
Q this code
10
0.68833
0.68867
-0.10416
-0.10423
100
0.67698
0.67695
-0.10152
-0.10094
1000
0.44799
0.44848
0.04992
0.04990
2000
0.29301
0.29311
0.0089989
0.00887
(nm)
*K. F. Evans and G. L. Stephens, “A new polarized atmospheric radiative transfer model,” J. Quant.
Spectrosc. Radiat Transfer. 46, 413-423, (1991).
II Escuela de Optica Biomedica, Puebla, 2011
Mueller Matrix -microspheres
solutions
x 10 -7
1
-0.2
0.5
-0.2
0
0.2
-0.2
0
0
0.2
0
-0.5
0.2
-1 -0.2
0
0
0.2
x 10 -7
1
-0.2
0.5
-0.2
0
0.2
-0.2
x 10 -7
1
-0.2
0.5
0
0
0.2
0
0
0.2
-7
-0.2
0
0.2
-0.2
0
0.2
0
0.2
-0.2
0
0
0
0
0.2
0.2
0
0.2
-0.5
0.2
-1 -0.2
0
0
0
0
0.2
0
0.2
0
-0.5
0
0.2
0
0.2
0.5
0
-0.5
0
0.2
0.5 10
0.5
0 20
0
20
40
-7
0
0.2
Monte Carlo
50
100
150
200
150
200
150
200
150
200
50 100 150 200
0.2
50 100 150 200
50
50
50
50
100
150
200
100
150
200
100
150
200
50 100 150 200
50 100 150 200
50 100 150 200
50 100 150 200
50
100
50
100
50
100
50
100
150
150
150
150
200
200
200
200
50 100 150 200
-7
50 100 150 200
50 100 150 200
50 100 150 200
50
50
50
50
0
100
100
100
100
-0.5
150
200
150
200
150
200
150
200
0.5
0
50 100 150 200
-1
x 10
1
0
50 100 150 200
100
150
200
-7
-0.5
x 10
1
-0.2
0.5
-0.5
0.2
-1 -0.2
50
100
-1
x 10
1
0
50
100
-1
-7
30
-0.5
40
-1
50
100
x 10 -7
1
0
-0.5
0.2
-1 -0.2
x 10
1
0
-0.5
0.2
-1 -0.2
0.5
0
-0.5
0.2
-1 -0.2
0
-7
x 10
1
-0.2
0.5
x 10 -7
1
x 10 -7
1
-0.2
0.5
0
x 10
1
-0.2
0.5
-0.5
0.2
-1 -0.2
0
-7
x 10
1
-0.2
0.5
-0.5
0.2
-1 -0.2
0
0
-0.5
0.2
-1 -0.2
x 10
1
-0.2
0.5
-7
-0.2
0
-0.5
0.2
-1 -0.2
x 10 -7
1
-0.2
0.5
0
-0.5
0.2
-1 -0.2
x 10 -7
1
-0.2
0.5
-1
50 100 150 200
50 100 150 200
50 100 150 200
50 100 150 200
Experimental*
II Escuela de Optica Biomedica, Puebla, 2011
*Cameron et al.
Mueller Matrix-asymmetric
illumination
M icrosphere solution
q
z
Polarizer
y
x
Polarizer
Laser
Camera
II Escuela de Optica Biomedica, Puebla, 2011
Mueller Matrix - microspheres
solutions
m11
m12
m11
m13
m12
m13
1
1
1
1
1
1
2
2
2
2
2
2
3
3
3
3
3
3
4
1
2
3
4
4
1
m21
2
3
4
4
1
m22
2
3
4
4
2
m21
m23
4
4
2
m22
4
4
1
1
1
1
1
1
2
2
2
2
2
2
3
3
3
3
3
3
4
4
4
1
2
3
4
1
m31
2
3
4
1
m32
2
3
4
4
2
m31
m33
4
4
2
m32
4
4
1
1
1
1
1
2
2
2
2
2
2
3
3
3
3
3
3
4
4
4
2
3
4
1
2
3
4
Monte Carlo
1
2
3
4
4
2
4
4
2
4
4
Experimental
II Escuela de Optica Biomedica, Puebla, 2011
4
2
4
m33
1
1
2
m23
2
4
Side View Experiments
0.5 µm
microspheres
solution
200 ml glass
container
7x7x4 cm
7 cm
He-Ne
Laser 633nm
1cm
Polarizer
0.25cm
Analyzer
12 bit
digital camera
G4
computer
II Escuela de Optica Biomedica, Puebla, 2011
Results -1
Parallel
experimental
Parallel
Monte Carlo
II Escuela de Optica Biomedica, Puebla, 2011
Results -2
Parallel
Max Parallel
II Escuela de Optica Biomedica, Puebla, 2011
Tomorrow
• Euler and Quaternion based Monte Carlo
• How to run the codes
• Application
II Escuela de Optica Biomedica, Puebla, 2011

similar documents