Large-Scale Inverse Problems in Imaging: Two Case Studies Open Access

Knepper, Sarah (2011)

Permanent URL: https://etd.library.emory.edu/concern/etds/4b29b661z?locale=en
Published

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

Rights statement
  • Permission granted by the author to include this thesis or dissertation in this repository. All rights reserved by the author. Please contact the author for information regarding the reproduction and use of this thesis or dissertation.
School
Department
Degree
Submission
Language
  • English
Research Field
Keyword
Committee Chair / Thesis Advisor
Committee Members
Last modified

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