Transcript New Approach to Solving the Radiative Transport Equation

New Approach to Solving the
Radiative Transport Equation
and Its Applications
George Y. Panasyuk
Bioengineering
UPenn, Philadelphia
[email protected]
Outline of the talk
• New approach for solving the RTE in a 3D macroscopically
homogeneous medium;
• Application of the method to:
(1) Calculation of the RTE Green’s function for the case
of free boundary;
(2) Generation of forward data for an inverse problem in
optical tomography.
Spectral Method
( zI  V) x  b
M different values of parameter Z; V is a N  N
Solve system of N equations:
for
| x,| b are vectors of length N
“Naïve” approach:
zI+V  W ( z)
and solve
W( z) x  b
for each M values of Z with
Computational complexity: M  N 3
matrix
Spectral method:
1) Find eigenvectors
eigenvalues
2) For every Z,
| n and
vn of V
| n n | b
| x  
n
z  vn
Computational complexity:
N3  M  N2
M  N3
Spectral Method for the RTE
2
ˆ
ˆ
ˆ
ˆ
ˆ


RTE: (s   t ) I (r, s)   s  P(s, s )I (r, s )d sˆ   (r, sˆ )
How can we write this equation in the form ( z  W ) I   ?
We can try to expand I (r,sˆ) into a 3D Fourier integral and into the
basis of spherical harmonics Ylm ( ,  )  Ylm (sˆ) the lab frame:
I (r, sˆ)    eikr d3kIlm (k)Ylm (sˆ) 
l ,m
ˆ
ˆ
k R|I(k)
 +S|I(k)
 |  (k),
(x) ˆ (y) ˆ (z)
ˆ
ˆ
R=(R , R , R )
( x)
*
Rlm

sin

cos

Y
,l m
lm ( ,  )Yl m ( ,  ) sin  d d ,

etc.,
S is diagonal
Rotated Reference Frames
To avoid k-dependence, use spherical harmonics Y (sˆ, kˆ ) defined in a reference
frame whose z-axis is aligned with the direction of k (“rotated” frames):
z

x
y'
k
z'
y
Ylm (sˆ, kˆ ) 
l

m l
Dml m (k ,k , 0)Ylm (sˆ)
Wigner D-functions
Euler angles = polar angles
 
x'
I (r, sˆ)    Ilm (k)Ylm (sˆ, kˆ )eikr d3k
l ,m
ikAˆ I (k)  Sˆ I (k)   (k)
Slm,lm  llmm Sl
- diagonal, Sl  a  s (1  Pl ), Pl  g l In HG model
k | k| is the spectral parameter,
Alm,lm   mm Dll
where D  is
ll
a tridiagonal real symmetric matrix with eigenfunctions Ψn and eigenvalues λn
1
I lm (k ) 
Sl

n
 (r , sˆ)   (r  r0 ) (sˆ  sˆ0 )  I 
 I (r , sˆ; r0 , sˆ0 ) 

 
m l ,l |m|
lm  n  n 
1  ik n
Green’s function of RTE
Ylm (sˆ, sˆ0 ) llm (r , r0 )Yl*m ( Rˆ , sˆ0 ), R  r  r0
Analytical dependences on all variables
Details: J.Phys.A 39, 115 (2006)
Evanescent Waves and the BVP
I kˆ,M ,n (r , sˆ)  e
 kˆ r / n
ˆ) lM  / S1/ 2 ,
ˆ
Y
(
s
,
k
 lM
n
l
kˆ  kˆ  1
l

Evanescent waves: kˆ (q)  inq  zˆnQ(q), q  zˆ  0,
I q,M ,n (r , sˆ)  eiq
Q(q ) z
Q(q)  q 2  1/ n2
ˆ ) lM  / S 1/ 2
ˆ
Y
(
s
,
k
 lM
n
l
l
vacuum
Z=0
ρ
Half-space z > 0
medium
z
I (  , z  0; sˆ)  Iinc (  , sˆ),
sˆ  zˆ  0
Solution of the half-space BVP:
2

ˆ
I (r , s)   d q Fq ,M ,n I q ,M ,n (r , sˆ)
n, M
J.Phys.A 39, 115 (2006)
Point uni-directional (sharply-peaked) source in an infinite
zˆ
medium placed at r0 = (0, 0, 0) and illuminating in the direction.
Forward and backword propagation: r = (0, 0, z),   zˆ, sˆ
zˆ, sˆ
zˆ
 
Angular dependence
ofzˆ, sˆthe specific intensity for forward (a) and backward (b)
propagation obtained at lmax = 21, g = 0.98 and a /s = 6 · 10-5 . The distance to the
source z is assumed to be positive for forward propagation and negative
for backward propagation.
*
is the transport mean free path.
Point uni-directional (sharply-peaked) source in an infinite medium
placed at r0 = (0, 0, 0) and illuminating in the zˆ -direction.
Off-axis case: r = (0, y, 0)
S0=
Two cases:
a) s is in the yz plane
z
z
  zˆ, sˆ
s

b) s is in the xy plane:
y
s

x
  yˆ , sˆ
*
*
*
*
*
*
Angular distribution of specific intensity for off-axis propagation (small and average absorption)
Parameters: g = 0.98 and a /s = 6 · 10-5 (a), (b), a /s = 0.03 (c), (d).
3.0
 2
(l ) I 10
16
10-16
lmax = 1
(l  )2 I
lmax = 3
lmax = 34
lmax =10
lmax = 34
0
10-18
-3.0
10-20
0
π
α, rad
2π
0
α0
π
α, rad 2π
Convergence of the specific intensity with lmax (left) and the converged
result at lmax = 34 (right); r = (0, 26l*, 0), φ = 0, g = 0.98, μa/μs = 0.2
s0
0 y
(a) Dependence of the position of maximum  0 on the distance to the source, y, for
physiological parameters: g  0.98 and a /s = 6 · 10-5 .
(b) Schematic illustration of typical "photon trajectories" that correspond to maxima in
specific intensity.
Application to optical tomography
-Recover absorption coefficient  (r )  0   (r ) of
inhomogeneous medium from multiple measurements with different
source-detector pairs, encoded in data function  (r1 , r2 ) :
1

 (r1 , r2 )   d r G(r1 , r )G(r , r2 ) 
3

3
*2

r G(r1 , r ) r G(r , r2 )  (r )

-G is Green’s function within DA
-l* is the transport free path
CCD
-Based on DA when lowest order correction in l* are
taken into account,
-Applicable when diffusion theory breaks down
(thin samples, near boundaries or sources, etc)
APL, 87, 101111 (2005)
1.0
corr.
no corr.
δα(x)/ δα(0)
slab thickness = 0.5cm
0.5
2Δx
Δx = 0.2l*
0
-3.0
-1.5
0
1.5
x/l*
3.0
1D profiles of reconstructed absorption coefficient α(x) = α0 + δα(x)
of a point absorber using corrected (red) and uncorrected (green) DA
Data function  (r1, r2 ) was simulated by the MRRF for the RTE
APL, 87, 101111 (2005)
CONCLUSIONS
●
●
●
The method of rotated reference frames takes advantage
of all symmetries of the RTE (symmetry with respect to
rotations and inversions of the reference frame).
The angular and spatial dependence of the solutions is
expressed in terms of analytical functions.
The analytical part of the solution is of considerable
mathematical complexity. This is traded for relative
simplicity of the numerical part. We believe that we have
reduced the numerical part of the computations to the
absolute minimum allowed by the mathematical structure
of the RTE
Co-Authors:
Vadim A. Markel
John C. Schotland
Publications:
1. V.A.Markel, "Modified spherical harmonics method for solving the radiative
transport equation," Letter to the Editor, Waves in Random Media 14(1), L13-L19
(2004).
2. G.Y.Panasyuk, J.C.Schotland, and V.A.Markel, "Radiative transport equation in
rotated reference frames," Journal of Physics A, 39(1), 115-137 (2006).
3. G.Y.Panasyuk, V.A.Markel, and J.C.Schotland, Applied Physics Letters 87, 101111
(2005).
Available on the web at
http://whale.seas.upenn.edu/vmarkel/papers.html