Numerical Modeling for Inverse Problems
Download
Report
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 ) Ji ...
• Rewriting: i y f (i ) y f (i 1 ) Ji i 1 Ji
• Substituting:
• Update equation:
J T i 0
J T ( i1 J ) 0
[J T J]i J T i1
[J T J I ]i J T i1
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!