Large-Scale Inverse Problems in Imaging: Two Case Studies Pubblico
Knepper, Sarah (2011)
Abstract
Solving inverse problems is an important part of scientific
computing. As computers be-
come more powerful, solutions to increasingly larger problems are
sought, allowing for more
accurate representations of real-world applications. We consider
solving large-scale inverse
problems, ranging from linear to fully nonlinear. We look at
aspects common to inverse
problems, such as their ill-posedness, and see how regularization
can help produce mean-
ingful results. We discuss a number of different methods for
solving while providing regu-
larization. One such technique is to solve using an iterative
method but stop the iterations
early, before convergence is fully achieved.
Iterative solvers are particularly useful for large-scale inverse
problems as computations
can be done in parallel. Trilinos is a mathematical software
library for solving problems
coming from many fields of scientific computing. One particular
package, Belos, provides
both an abstract framework and concrete implementations of various
iterative solvers. We
have implemented two additional solvers within the Belos framework,
LSQR and MRNSD,
which can be used to solve linear inverse problems.
We then consider two different case studies, where we wish to solve
a large-scale linear
inverse problem. In the first study, we want to remove patient
motion blur from positron
emission tomography (PET) images when motion information is tracked
and recorded dur-
ing the scan. We describe how this problem can be formulated as a
linear equation, then
we solve it using the solvers we implemented. We also look at a
number of results, seeing
how the reconstruction improves as more motion information is
included in our model.
The second case study comes from the field of adaptive optics. Here
we wish to determine
the distortion caused by the atmosphere when imaging using
ground-based telescopes. Sen-
sors are able to obtain noisy estimates of the gradients of the
distortion, resulting in a Kro-
necker product-structured linear least squares problem. We describe
a solving method that
employs Tikhonov-type regularization by exploiting properties of
the Kronecker product
and utilizing the generalized singular value decomposition (GSVD).
Our approach includes
constructing a preconditioner off-line and then applying a few
iterations of preconditioned
LSQR.
Table of Contents
1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . 2
1.1.1 Model Problems . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . 2
1.2 Mathematical Modelling and Analysis . . . . . . . . . . . . . .
. . . . . . . 4
1.2.1 Linear Inverse Problems . . . . . . . . . . . . . . . . . . .
. . . . . . 4
1.2.2 Separable Nonlinear Inverse Problems . . . . . . . . . . . .
. . . . . 16
1.2.3 Nonlinear Inverse Problems . . . . . . . . . . . . . . . . .
. . . . . . 20
1.3 Outline of Work . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 23
1.4 Contributions . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . 23
2 Large-Scale Software for Inverse Problems 25
2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . 25
2.2 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 25
2.3 Overview of Trilinos . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . 26
2.3.1 Petra Model . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . 27
2.3.2 Teuchos Toolkit . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . 28
2.3.3 Belos Framework . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . 28
2.4 LSQR . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . 29
2.4.1 Implementation Details of LSQR . . . . . . . . . . . . . . .
. . . . . 33
2.5 MRNSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 33
2.5.1 Implementation Details of MRNSD . . . . . . . . . . . . . . .
. . . . 34
2.6 Least Error Convergence Test . . . . . . . . . . . . . . . . .
. . . . . . . . . 35
2.7 Remarks and Future Directions . . . . . . . . . . . . . . . . .
. . . . . . . . 36
3 Case Study 1: Positron Emission Tomography Application 37
3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . 37
3.2 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 38
3.3 Methodology of Our Approach . . . . . . . . . . . . . . . . . .
. . . . . . . 39
3.3.1 Motion Detection . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . 39
3.3.2 Construction of the Matrix . . . . . . . . . . . . . . . . .
. . . . . . 40
3.3.3 Iterative Deblurring . . . . . . . . . . . . . . . . . . . .
. . . . . . . 42
3.4 Implementation Details . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 43
3.4.1 Memory Requirements . . . . . . . . . . . . . . . . . . . . .
. . . . . 43
3.4.2 Scalability Analysis . . . . . . . . . . . . . . . . . . . .
. . . . . . . 46
3.4.3 Testbed . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . 49
3.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . 49
3.5.1 Effect of Scalar Precision . . . . . . . . . . . . . . . . .
. . . . . . . 50
3.5.2 Effect of Interval Choice . . . . . . . . . . . . . . . . . .
. . . . . . . 53
3.5.3 Effect of Patient Motion . . . . . . . . . . . . . . . . . .
. . . . . . . 55
3.6 Remarks and Future Directions . . . . . . . . . . . . . . . . .
. . . . . . . . 58
4 Case Study 2: Adaptive Optics Application 61
4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . 61
4.2 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 63
4.3 Methodology of Our Approach . . . . . . . . . . . . . . . . . .
. . . . . . . 63
4.3.1 Mathematical Background . . . . . . . . . . . . . . . . . . .
. . . . . 64
4.3.2 Wavefront Reconstruction Using TSVD-Type Regularization . . .
. 66
4.3.3 Wavefront Reconstruction Using Tikhonov-Type Regularization .
. . 68
4.4 Implementation Details . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 73
4.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . 73
4.5.1 Effect of Differing alpha Values . . . . . . . . . . . . . .
. . . . . . . . . 75
4.5.2 Using Square-Aperture Preconditioner on Masked Problems . . .
. . 75
4.6 Remarks and Future Directions . . . . . . . . . . . . . . . . .
. . . . . . . . 79
5 Concluding Remarks 80
A Trilinos Code 82
A.1 Code for LSQR . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 82
A.1.1 LSQRSolMgr.hpp Code . . . . . . . . . . . . . . . . . . . . .
. . . . 82
A.1.2 LSQRIter.hpp Code . . . . . . . . . . . . . . . . . . . . . .
. . . . . 96
A.1.3 LSQRStatusTest.hpp Code . . . . . . . . . . . . . . . . . . .
. . . . 107
A.2 Code for MRNSD . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 112
A.2.1 MRNSDSolMgr.hpp Code . . . . . . . . . . . . . . . . . . . .
. . . . 112
A.2.2 MRNSDIter.hpp Code . . . . . . . . . . . . . . . . . . . . .
. . . . . 124
A.2.3 MRNSDStatusTest.hpp Code . . . . . . . . . . . . . . . . . .
. . . . 133
A.3 Code for Least Error Status Test . . . . . . . . . . . . . . .
. . . . . . . . . 137
A.3.1 LeastErrorStatusTest.hpp Code . . . . . . . . . . . . . . . .
. . . . . 137
A.4 Code for PET Application . . . . . . . . . . . . . . . . . . .
. . . . . . . . . 141
A.4.1 HRRT.hpp Code . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . 141
A.4.2 HRRTmain.cpp Code . . . . . . . . . . . . . . . . . . . . . .
. . . . 160
A.4.3 Example XML File . . . . . . . . . . . . . . . . . . . . . .
. . . . . . 160
A.5 Code for AO Application . . . . . . . . . . . . . . . . . . . .
. . . . . . . . 161
A.5.1 AOOperator.hpp . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . 161
A.5.2 AOPreconditioner.hpp . . . . . . . . . . . . . . . . . . . .
. . . . . . 166
A.5.3 AOmain.hpp . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . 169
A.5.4 AOmain.cpp . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . 176
A.5.5 Example XML File . . . . . . . . . . . . . . . . . . . . . .
. . . . . . 176
About this Dissertation
School | |
---|---|
Department | |
Degree | |
Submission | |
Language |
|
Research Field | |
Parola chiave | |
Committee Chair / Thesis Advisor | |
Committee Members |
Primary PDF
Thumbnail | Title | Date Uploaded | Actions |
---|---|---|---|
Large-Scale Inverse Problems in Imaging: Two Case Studies () | 2018-08-28 13:48:09 -0400 |
|
Supplemental Files
Thumbnail | Title | Date Uploaded | Actions |
---|---|---|---|
SingularValuesDeriv2B10.eps () | 2018-08-28 13:54:43 -0400 |
|
|
reductionNNMRNSD.eps () | 2018-08-28 13:56:05 -0400 |
|
|
FullDoubleErrorIntervals.eps () | 2018-08-28 13:57:25 -0400 |
|
|
RelDifferenceSingularValuesDeriv2B20.eps () | 2018-08-28 13:58:52 -0400 |
|
|
BigSigma.eps () | 2018-08-28 14:00:13 -0400 |
|
|
AOPreconditioner.hpp () | 2018-08-28 14:02:25 -0400 |
|
|
reductionNNLSQR.eps () | 2018-08-28 14:03:41 -0400 |
|
|
MRNSDSP.eps () | 2018-08-28 14:05:06 -0400 |
|
|
TikhonovRecon256.eps () | 2018-08-28 14:06:26 -0400 |
|
|
ClassDiagram2.pdf () | 2018-08-28 14:07:37 -0400 |
|
|
ItersMRNSD.eps () | 2018-08-28 14:09:18 -0400 |
|
|
HRRTmain.cpp () | 2018-08-28 14:10:48 -0400 |
|
|
AOmain.cpp () | 2018-08-28 14:12:46 -0400 |
|
|
reductionTriMRNSD.eps () | 2018-08-28 14:13:45 -0400 |
|
|
timings560mini.eps () | 2018-08-28 14:15:14 -0400 |
|
|
SPDPIterationsLSQR.eps () | 2018-08-28 14:16:46 -0400 |
|
|
SingularValuesPlotN64.eps () | 2018-08-28 14:18:18 -0400 |
|
|
MRNSDDP.eps () | 2018-08-28 14:19:46 -0400 |
|
|
MotionInformation7Intervals.eps () | 2018-08-28 14:22:27 -0400 |
|
|
NoisySy256.eps () | 2018-08-28 14:23:12 -0400 |
|
|
SPDPIterationsMRNSD.eps () | 2018-08-28 14:23:49 -0400 |
|
|
BilinearInterp.eps () | 2018-08-28 14:26:15 -0400 |
|
|
NearestNeighborInterp.eps () | 2018-08-28 14:27:53 -0400 |
|
|
PET.xml () | 2018-08-28 14:30:04 -0400 |
|
|
ItersTolLevelPrec.eps () | 2018-08-28 14:31:25 -0400 |
|
|
timings2500mini.eps () | 2018-08-28 14:33:01 -0400 |
|
|
RelErrLowMotion.eps () | 2018-08-28 14:34:35 -0400 |
|
|
NoisySx256.eps () | 2018-08-28 14:35:39 -0400 |
|
|
MaskedRelErrTolLevelPrec.eps () | 2018-08-28 14:36:49 -0400 |
|
|
MotionMid.eps () | 2018-08-28 14:37:57 -0400 |
|
|
RelErrLSQR.eps () | 2018-08-28 14:39:23 -0400 |
|
|
SingularValuesDeriv2B20.eps () | 2018-08-28 14:40:50 -0400 |
|
|
FullDifferenceIterationIntervals.eps () | 2018-08-28 14:42:13 -0400 |
|
|
LeastErrorStatusTest.hpp () | 2018-08-28 14:45:06 -0400 |
|
|
RelDifferenceSingularValuesDeriv2B50.eps () | 2018-08-28 14:46:09 -0400 |
|
|
MotionAll.eps () | 2018-08-28 14:47:41 -0400 |
|
|
FullFloatErrorIntervals.eps () | 2018-08-28 14:49:16 -0400 |
|
|
RelErrIdealEqualTriMRNSD.eps () | 2018-08-28 14:50:46 -0400 |
|
|
TSVDRecon256.eps () | 2018-08-28 14:52:25 -0400 |
|
|
RelErrHighMotion.eps () | 2018-08-28 14:53:53 -0400 |
|
|
LSQRIter.hpp () | 2018-08-28 14:55:21 -0400 |
|
|
RelSpreadSingularValuesDeriv2.eps () | 2018-08-28 14:56:48 -0400 |
|
|
RelErrIdealEqualNNLSQR.eps () | 2018-08-28 14:58:17 -0400 |
|
|
LSQRDP.eps () | 2018-08-28 15:00:10 -0400 |
|
|
LSQRStatusTest.hpp () | 2018-08-28 15:01:38 -0400 |
|
|
RelDifferenceSingularValuesDeriv2B10.eps () | 2018-08-28 15:03:11 -0400 |
|
|
RelErrMRNSD.eps () | 2018-08-28 15:05:24 -0400 |
|
|
AOOperator.hpp () | 2018-08-28 15:06:47 -0400 |
|
|
ItersLowMotion.eps () | 2018-08-28 15:08:16 -0400 |
|
|
SingularValuesDeriv2.eps () | 2018-08-28 15:10:30 -0400 |
|
|
LSQRSP.eps () | 2018-08-28 15:12:43 -0400 |
|
|
MotionHigh.eps () | 2018-08-28 15:14:09 -0400 |
|
|
ItersMidMotion.eps () | 2018-08-28 15:15:23 -0400 |
|
|
RelErrIdealEqualTriLSQR.eps () | 2018-08-28 15:16:37 -0400 |
|
|
MaskedRelErrTolLevelUnprec.eps () | 2018-08-28 15:19:26 -0400 |
|
|
ItersTolLevelUnprec.eps () | 2018-08-28 15:21:46 -0400 |
|
|
MotionLow.eps () | 2018-08-28 15:23:25 -0400 |
|
|
timings2500small.eps () | 2018-08-28 15:24:55 -0400 |
|
|
FullFloatIterationIntervals.eps () | 2018-08-28 15:26:18 -0400 |
|
|
RelErrMidMotion.eps () | 2018-08-28 15:27:44 -0400 |
|
|
SPDPRelErrLSQR.eps () | 2018-08-28 15:29:11 -0400 |
|
|
MRNSDStatusTest.hpp () | 2018-08-28 15:30:41 -0400 |
|
|
ItersHighMotion.eps () | 2018-08-28 15:32:00 -0400 |
|
|
BlurredLargeHeadSlice24PatN.eps () | 2018-08-28 15:34:17 -0400 |
|
|
AnnularMask.eps () | 2018-08-28 15:37:12 -0400 |
|
|
RelErrIdealEqualNNMRNSD.eps () | 2018-08-28 15:38:00 -0400 |
|
|
timings560small.eps () | 2018-08-28 15:40:50 -0400 |
|
|
CircularMask.eps () | 2018-08-28 15:42:08 -0400 |
|
|
MRNSDIter.hpp () | 2018-08-28 15:43:50 -0400 |
|
|
AO.xml () | 2018-08-28 15:45:20 -0400 |
|
|
MRNSDSolMgr.hpp () | 2018-08-28 15:46:47 -0400 |
|
|
AOmain.hpp () | 2018-08-28 15:48:15 -0400 |
|
|
AlphaIterationsN256Bold.eps () | 2018-08-28 15:49:42 -0400 |
|
|
reductionTriLSQR.eps () | 2018-08-28 15:51:11 -0400 |
|
|
SPDPRelErrMRNSD.eps () | 2018-08-28 15:53:50 -0400 |
|
|
HRRT.hpp () | 2018-08-28 15:54:32 -0400 |
|
|
SingularValuesDeriv2B50.eps () | 2018-08-28 15:57:15 -0400 |
|
|
FullDoubleIterationIntervals.eps () | 2018-08-28 15:58:31 -0400 |
|
|
LSQRSolMgr.hpp () | 2018-08-28 16:00:19 -0400 |
|
|
MotionInformation14Intervals.eps () | 2018-08-28 16:01:50 -0400 |
|
|
FullDifferenceErrorIntervals.eps () | 2018-08-28 16:03:26 -0400 |
|
|
IterationsCountN256Bold.eps () | 2018-08-28 16:05:10 -0400 |
|
|
ItersLSQR.eps () | 2018-08-28 16:06:45 -0400 |
|