Optical Tomography Based on the Method of Rotated Reference Frames - - PowerPoint PPT Presentation

Optical Tomography Based on the Method of Rotated Reference Frames

Manabu Machida (Dept. of Bioengineering, Univ. of Pennsylvania)

[Applied Inverse Problems, Vienna, July 20, 2009]

In collaboration with George Y. Panasyuk (U. Penn), Zheng-Min Wang (U. Penn), Vadim A. Markel (U. Penn), and John C. Schotland (U. Penn)

Outline

• Introduction
• Green’s function
• Inverse problem
• Simulation
• Experiment

Optical Tomography

e.g.,

3D slab geometry

• s,z = 0

( )

• d,z = L

( )

z 104 sources 105 detectors hs hd

Transport (RTE) Regime

Maxwell eqs Transport eq Diffusion approx

• Multiple scattering
• Low scattering
• Near sources
• Near boundaries
• Thin samples
• etc.

depth

OCT DOT

Needs to use the radiative transport equation (RTE)

☛ No correlation among photons ☛ No time dependence ☛ Background is homogeneous

O. Dorn, Inverse Problems (1998) A. D. Klose and A. H. Hielscher, Med. Phys. (1999) S. Wright, M. Schweiger, and S. R. Arridge,

• Meas. Sci. Technol. (2007)

T. Tarvainen, M. Vauhkonen, and S. R. Arridge,

• J. Quant. Spec. Rad. Trans. (2008)

P. Gonzalez-Rodriguez and A. D. Kim, Inverse Problems (2009)

OT by Radiative Transport Equation

RTE in the 3D slab geometry

Outline

• Introduction
• Green’s function
• Inverse problem
• Simulation
• Experiment

[ : scattering and absorption]

Green’s Function for the RTE

ˆ LI0 ,z,ˆ s

( ) =

d2 s A ˆ s,ˆ s

( )I0 ,z,ˆ

s

( )

• µt = µa + µs

A ˆ s ˆ s

( )

I0 ,z,ˆ s

( )

µs, µa

: specific intensity : phase function

ˆ s + µt µs ˆ L

• I0 ,z,ˆ

s

( ) =

( ) z

( ) ˆ

s ˆ z

( )

µt and µs = const

g = d2 s ˆ s ˆ s

( )A ˆ

s,ˆ s

( )

• : scattering asymmetry parameter

= G ,z,ˆ s;

0,z0,ˆ

s0

( )

= x,y

( )

Boundary Conditions

I0 ,zb,ˆ

s, ˆ s

( ) = Rn

ˆ n ˆ s

( )I0 ,zb, ˆ

s, ˆ s

( ),

ˆ n ˆ s < 0 Rn x

( ) =

r

p 2 + r s 2

2 x xc 1 x < xc

• x0 =

1 n2 1 x2

( )

xc = n2 1 n r

p = x x0

x + x0 r

s = x0 x

x0 + x ˆ n ˆ n ˆ s

Snell and Fresnel

ˆ s = ˆ

s, ˆ s

( )

1 n 1

V. A. Markel, Waves Random Media 14 (2004) L13 G. Panasyuk, J. C. Schotland, and V. A. Markel,

• J. Phys. A: Math. Gen. 39 (2006) 115

J. C. Schotland and V. A. Markel,

• Inv. Prob. Imag. 1 (2007) 181

M. M., G. Y. Panasyuk, J. C. Schotland, and V. A. Markel, submitted.

Method of Rotated Reference Frames

Green’s function for the 3D slab geometry

ˆ s G ,z,ˆ s;

0,0, ˆ

z

( ) + µtG ,z,ˆ

s;

0,0, ˆ

z

( ) =

µs d2 s A ˆ s,ˆ s

( )

• G ,z,ˆ

s ;

0,0, ˆ

z

( ) + ( ) z

( ) ˆ

s ˆ z

( )

Structure of the Green’s Function

G ,z,ˆ s;

0,0, ˆ

z

( ) =

d2q 2

( )

2 F qµ +

( )

( )Iqµ

+

( ) ,z,ˆ

s

( ) + F

qµ

• ( )

( )Iqµ

• ( ) ,z,ˆ

s

( )

• µ
• G ,0,ˆ

s

( ) =

( ) ˆ

s ˆ z

( ) + Rn ˆ

s ˆ z

( )G ,0, ˆ

s, ˆ s

( )

ˆ s ˆ z > 0

G ,L,ˆ s

( ) = Rn ˆ

s ˆ z

( )G ,L, ˆ

s, ˆ s

( )

ˆ s ˆ z < 0

MRRF Boundary Conditions

◆ Numerical diagonalization of a tri-diagonal matrix. ◆ Dependences on , , and are analytical.

Basis Modes

z ˆ s

exp iq Qµ q

( )z

• Ylm ˆ

s; ˆ k

( ) =

Yl

• m

ˆ s

( )D

m m l

ˆ

k, ˆ k,0

( )

• m =l

l

• Ylm ˆ

s

( ) Iqµ

+

( ) ,z,ˆ

s

( ), Iqµ

• ( ) ,z,ˆ

s

( )

Spherical harmonics in a rotated reference frame Plane-wave decomposition

• k = iq ± Qµ q

( ) ˆ

z

ˆ s Iqµ

±

( ) ,z,ˆ

s

( ) + µtIqµ

±

( ) ,z,ˆ

s

( ) µs d2

s A ˆ s,ˆ s

( )

• Iqµ

±

( ) ,z,ˆ

s

( ) = 0

Green’s Function

G ,z,ˆ s;

0,0, ˆ

z

( ) =

d2q 2

( )

2 e iq 0

( )

• Ylm ˆ

s

( )e

imqklm q,z;n

( )

m=l l

• l=0
• klm q,z;n

( ) =

l n M

( )

l

Mn

• dmM

l

i qn M

( )

( )

• e

QMn q

( )z fMn

+

( ) q;n

( ) + e

QMn q

( ) Lz ( ) 1

( )

l+m+ M fMn

• ( ) q;n

( )

• l = µa + µs 1 gl

( )

QMn q

( ) =

q2 + µ

2

[Henyey-Greenstein phase function]

Green’s Function

G ,z,ˆ s;

0,0, ˆ

z

( ) =

d2q 2

( )

2 e iq 0

( )

• Ylm ˆ

s

( )e

imqklm q,z;n

( )

m=l l

• l=0
• klm q,z;n

( ) =

l n M

( )

l

Mn

• dmM

l

i qn M

( )

( )

• e

QMn q

( )z fMn

+

( ) q;n

( ) + e

QMn q

( ) Lz ( ) 1

( )

l+m+ M fMn

• ( ) q;n

( )

• l = µa + µs 1 gl

( )

QMn q

( ) =

q2 + µ

2

[Henyey-Greenstein phase function] Wigner’s d-function eigenvector eigenvalue b.c. spherical harmonics

Outline

• Introduction
• Green’s function
• Inverse problem
• Simulation
• Experiment

Structure of Our Method

µa(,z) = µa + µa(,z) I0 I G

• d,L, ˆ

z;,z,ˆ s

( )µa ,z

( )G ,z,ˆ

s;

s,0, ˆ

z

( )

µa ,z

( )

Inverse Problem Forward Problem

Signal Background

I0

d,L,ˆ

s;

s

( ) µa

I

d,L,ˆ

s;

s

( ) µa ,z

( )

Absorbing inhomogeneity

Markel and Schotland,

• Phys. Rev. E (2004)

Schotland and Markel,

• Inv. Prob. Imag. (2007)

Radiative Transport Equation

ˆ s + µt + µa ,z

( ) µs ˆ

L

• I ,z,ˆ

s

( ) = ,z,ˆ

s

( )

ˆ LI ,z,ˆ s

( ) =

d2 s A ˆ s,ˆ s

( )I ,z,ˆ

s

( )

• ,z,ˆ

s

( )

s

( ) z

( ) ˆ

s ˆ z

( )

ˆ s + µt µs ˆ L

• I0 ,z,ˆ

s

( ) = ,z,ˆ

s

( )

Only and in the normal direction is detected through the aperture.

I

d,L, ˆ

z

( )

I0

d,L, ˆ

z

( )

source:

Born Approximation

d2dzd2s G

• d,L, ˆ

z;,z,ˆ s

( )µa ,z

( )G ,z,ˆ

s;

s,0, ˆ

z

( )

s, d

( ) I0

d,L, ˆ

z

( ) I

d,L, ˆ

z

( )

Data Function:

I

d,L, ˆ

z

( ) I0

d,L, ˆ

z

( )

• d2dzd2s G
• d,L, ˆ

z;,z,ˆ s

( )µa ,z

( )I ,z,ˆ

s

( )

I0 ,z,ˆ s

( )

Fourier Transform (1)

• qs,qd

( ) =

e

i qs s +qd d

( )

s, d

( )

d

• s
• µa q,z

( ) =

d2 eiqµa ,z

( )

• G ,z,ˆ

s;

0,0, ˆ

z

( ) =

d2q 2

( )

2 e iq 0

( )

• Ylm ˆ

s

( )e

imqklm q,z;n

( )

m=l l

• l=0

Fourier Transform (2)

• qs,qd

( )

1 hs

2hd 2

dz qs,qd,z

( )

µa qs + qd,z

( )

• qs,qd,z

( ) =

e

im qs +qd

( )klm qs,z;n

( )klm qd,L z;n ( )

lm

• qs = q

2 + p, qd = q 2 p q,p

( ) hs

2hd 2

q 2 + p, q 2 p

• K q,p,z

( ) p + q

2 ,p q 2 ,z

Fourier Transform (3)

q,p

( ) =

dz K q,p,z

( )

µa q,z

( )

• s,

d

( ) =

d2dzd2s G

• d,L, ˆ

z;,z,ˆ s

( )µa ,z

( )G ,z,ˆ

s;

s,0, ˆ

z

( )

= ˆ K µa

Singular Value Decomposition

• µa q,z

( ) =

E j,

( )

E j

2 j

• K * q,

p ,z

( )v j

• p

( )

• p
• v j

p

( )*

p

• q,p

( )

• µa = ˆ

K † ˆ K † = ˆ K * ˆ K ˆ K *

( )

1

M p

p q

( ) =

dz K q,p,z

( )K * q,

p ,z

( )

• M p

p q

( )v j

• p

( ) q

( )

• p
• = E j

2 q

( )v j

p

( ) q

( )

µa ,z

( ) =

d2q 2

( )

2 eiq

• µa q,z

( )

Regularizer Ill-posed

Outline

• Introduction
• Green’s function
• Inverse problem
• Simulation
• Experiment

Bars

L = 6*, hs = 0.2*, hd = 0.1*, µa = 0.05 cm-1, µs 1 g

( ) = 10 cm-1,

g = 0.9, * = 1 mm

[RTE] [DA]

x-y planes at depth z

Balls

z x

[RTE] [DA]

x-y planes at depth z

L = 6*, hs = 0.2*, hd = 0.1*, µa = 0.05 cm-1, µs 1 g

( ) = 10 cm-1,

g = 0.9, * = 1 mm

Outline

• Introduction
• Green’s function
• Inverse problem
• Simulation
• Experiment

Real World (1)

Laser CCD 3mm 3mm

Width of slab: L=10mm

Real World (2)

hs = 3.2 mm, hd = 0.36 mm, 292 sources 3972 detectors, µa = 0.03 cm-1, µs 1 g

( ) = 7 cm-1,

g = 0.65, * = 1.4 mm, n = 1.8

[RTE] [DA]

x-y planes at depth z

• Three-dimensional objects with arbitrary shapes

embedded in the slab can be reconstructed with the radiative transport equation.

• Computation time for reconstruction is 90 sec

with four-cpu parallel computation of 1.3 GHz Itanium-2.

• In experiments, to infrared sources are

collected by pixels of the CCD.

Summary

104 103 105

Diffusion Approximation

Q q

( ) =

q2 + 3µa *

absorbing boundary conditions

q,p

( ) =

dz KDA q,p,z

( )

µa q,z

( )

• KDA q,p,z

( ) DA p + q

2 ,p q 2 ,z

• DA q1,q2,z

( ) =

3 1 µa*

( )

4

• 2 sinh Q q1

( ) L z

( )

• sinh Q q2

( )z

• sinh Q q1

( )L

• sinh Q q2

( )L

Fitting

DA (experiment)

DA q1,q2,z

( ) =

3 µs 4

• 2

* +

( )

2 gb q1;0,z

( )gb q2;L,z ( )

gb q,z1,z2

( ) =

sinh Q L z2 z1

( )

• + Qcosh Q L z2 z1

( )

• sinh QL

[ ]+ 2Qcosh QL [ ]+ Q

( )

2 sinh QL

[ ]

Q q

( ) =

q2 + k2

GDA ,z;

0,0

( ) = 3

c* d2q 2

( )

2 gb q,z,0

( )e

iq 0

( )

• k = 0.02 / hd,

= 1.2hd

fitting

Rytov Approximation

s, d

( ) = G

d,L, ˆ

z;

s,0, ˆ

z

( )ln I0

d,L, ˆ

z

( )

I

d,L, ˆ

z

( )

s, d

( ) = I0

d,L, ˆ

z

( ) I

d,L, ˆ

z

( )

Data Function:

I

d,L, ˆ

z

( ) I0

d,L, ˆ

z

( )

exp d2dzd2s G

d,L, ˆ

z;,z,ˆ s

( )µa ,z

( )G ,z,ˆ

s;

s,0, ˆ

z

( )

• G

d,L, ˆ

z;

s,0, ˆ

z

( )

• cf., Born approximation: