Transcript Numerical Modeling for Inverse Problems

Numerical Modeling for
Image Reconstruction
Subha Srinivasan
11/2/2009
Definition of Inverse Problem..
• Definition: Given a distribution of sources
and a distribution of measurements at the
boundary dΩ, finding the tissue parameter
distribution within domain Ω.
• Expressed as x = F-1(y)
Ways of Solving Inverse Problems
• Back-projection methods
• Perturbation methods
• Non-linear optimization methods
Back-projection methods
patient
fan beam
detectors
x-rays
Filtered Back-projection method:
[measurements] = [attenuation op.] [object ]
[image] = [attenuation op.]T [filter] [measurements]
• Assumes that each
projection provides a
nearly independent
measurement of the
domain.
• Assumes that light
travels in a straight
line: not true with
tissue unless
scattering is isolated
Linear reconstruction for Change in
Optical Properties
• x = F-1(y) is a non-linear problem: can be linearized using
Taylor’s series expansion if initial estimate is close to actual
values:
y  f (0 )  f '(0 )(  0 )  ...
y  J
y  (y  y0 ), where y0  f (0 )
Jacobian
matrix
J
f ( 0 )

Reconstructing for changes rather than absolute values
Structure of Jacobian
C
J

  ln I 1  ln I 1

D1
D 2

 1
 1


D1
D 2
  ln I
 ln I 2
2

D1
D 2

J 
 2
 2

D1
D 2




  ln I NM  ln I NM
 D
D 2
1

 NM
  NM
 D1
D 2
 ln I 1
;
D NN
 1

;
D NN
 ln I 2

;
D NN
 2

;
D NN

Ф = I e-i(ωt+θ)
 = [a, D]
I = signal amplitude
θ = signal phase
Absorption coeff.
Diffusion Coeff.
 ln I 1
 a1
 1
 a1
 ln I 2
 a1
 2
 a1



 ln I NM
 ln I NM

;
D NN
 a1

 NM
;
D NN
 NM
 a1
 ln I 1
 a 2
 1
 a 2
 ln I 2
 a 2
 2
 a 2

 ln I NM
 a 2
 NM
 a 2
 ln I 1
 aNN
 1

 aNN
 ln I 2

 aNN
 2

 aNN





 ln I NM
 aNN
 NM
 aNN

















Shape of Jacobian
Calculated by:
1)Perturbation
Method
2)Direct Analytic
Jacobian
3) Adjoint method
Adjoint method for Jacobian Calculation
PDIRECT  PADJOINT
J a (s, j) 
J (s, j) 
s, j
a
s, j
.adj

 s  Adj
j
 s   Adj
j
i  Adj

  a     q Adj

c
Dehghani notes
Solving
• Linearizing change in intensity: born approximation
• Linearizing change in log intensity: Rytov
approximation
• Inverting J: large, under-determined and ill-posed:
some standard methods can be used
• Truncated SVD, Tikhonov regularization, Algebriac
reconstruction techniques (ART) & Conjugate
Gradient methods are commonly used
Terminology: Inverse Problem
• Ill-posed–Small changes in the data can cause large
changes in the parameters.
• Ill-conditioned–The condition number (ratio of
largest singular value to smallest singular value) is
large, which implies the inverse solution would not be
unique.
• Ill-determined–(or under-determined) The number
of independent equations are smaller than number of
unknowns.
Deriving Update Equation using
Least Squares Minimization
• Minimizing error functional:
2

 J T  0

• Setting derivative to zero:
• Taylor’s approximation
  y  f ( )
f (i )  f (i 1 )  Ji  ...
• Rewriting: i  y  f (i )  y  f (i 1 )  Ji  i 1  Ji
• Substituting:
• Update equation:
J T i  0
 J T ( i1  J )  0
[J T J]i  J T i1
[J T J   I ]i  J T i1
Assumptions of Levenberg-Marquardt
Minimization:
• JTJ is positive-definite
• Initial guess must be close to actual solution
• Update equation does not solve first-order
conditions unless α = 0
*Yalavarthy et. al., Medical Physics, 2007
Tikhonov Minimization:
Key idea is to introduce apriori assumptions about size and smoothness of desired solution:
L is dimensionless
common choice: L = I
(the identity matrix)
*Tikhonov et. al, 1977; Tarantola SIAM 2004.
*Yalavarthy et. al., Medical Physics, 2007
Tikhonov Minimization
Advantage:
• parameters within the minimization scheme
=> stability
Limitation:
• it requires a prior opinion about the noise
characteristics of the parameter and data
spaces (for λ)
Flow-chart for Iterative Image Reconstruction
Choosing Regularization: L-curve criterion
• Convenient graphical tool for displaying trade-off
between size of solution and its fit to the given data as
λ varies.
• λ can also be chosen empirically or based on
parameter/data values.
Hansen, ‘L-curve and its use in numerical treatment of inverse problems’
Reconstruction Results
• Simulated Measurements, 5% Noise
Recovery of Absorption
Recovery of Scattering
Spectral Image Reconstruction
nc
a ( )    i ( )ci
i 1
s '( )  A b
Data from Boulnois et al, Hale & Quarry,
figure from thesis Srinivasan et al
Spectral Image Reconstruction
c  (JspT Jsp  sp I)1 JspT sp
   1

   2
 sp  
 M
   k








2
 c1 
c 
 2
M 
c 
 cn 
A 
 
b 
 J c1, 1

 J
J sp   c1,  2
 M
 J c1,  k

J c2, 1 K
J cn, 1
J A, 1
J c2,  2 K
J cn,  2
J A,  2
M
M
J cn,  k
J A,  k
M
J c2,  k K
Relationships between Jsp & J can be obtained
Details, refer to Srinivasan et al, AO, 2005
J b, 1 

J b,  2 

M 
J b,  k 

Simulations show Reduced Cross-talk in spectral images
•
Data generated from a tumor-simulating phantom using FEM forward model,
with 1% random-Gaussian noise added.
HbT(μM)
StO2(%)
Water (%)
Scatt Ampl.
Scatt Power
True
Spectral
Conv.
• Spectral Method: Smoother Images; 15.3 % mean error compared to 43% (conv. Method).
• Reduced Cross-talk between HbO2 and water: from 30% (conv.) to 7% (spectral).
•Accuracy in StO2 accurate (<1% error)
Srinivasan et al, PhD thesis, 2005
Results from Image Reconstruction:
Experimental Data
Brooksby, Srinivasan et al, Opt Lett, 2005
References
• Gibson et al, Phy Med Bio: 50 : 2005: A review paper
• Paulsen et al, Med Phy: 22(6): 1995: first results from imagereconstruction in DOT
• Yalavarthy et al, Med Phy: 34(6): 2007: good explanation of math
• Brooksby et al, IEEE Journal of selected topics in quantum
electronics: 9(2): 2003: good reference for spatial priors
• Hansen: ‘Rank deficient and discrete ill-posed problems’:
SIAM: 1998: good reference for tikhonov/l-curve
• Srinivasan et al, Appl Optics: 44(10): 2005: reference for spectral
priors
• Press et al: ‘Numerical Recipes in Fortran 77’: II edition:
1992: great book for numerical folks!