Inverse problems, spring 2017
Question: What on earth is this picture about?
Consider the above photographic situation. You would like to have taken the left image: sharp and beautiful. However, due to misfocusing and high ISO setting in the camera, you ended up with the blurry and noisy picture on the right. Is there any chance to save the day? Well, one can apply deblurring, one of the topics of this course. Below you see three different deblurring methods with varying success. From left to right: basic regularization, total variation, and total generalized variation. (Thanks to Professor Kristian Bredies for sharing his amazing TGV algorithm!)
Results of the mid-course questionnaire (Feb 15). It seems that all is going pretty well in general. Let's discuss in the lecture about how to make finding material easier.
Lecture 16 (March 24)
Project work info. The goals of this session are
- Dividing students into project teams of two members each,
- Choosing the topic (tomography or deconvolution data, what is the focus of the project, what methods to try),
- Understanding the goals of Phase 1 and Phase 2 of the project,
- Agreeing on meeting times for reviewing Phase 1 of the projects.
Lecture 15 (March 22)
Real-data X-ray tomography. How to go from projection images to line integral data? What are the crucial properties of the variants of tomography, such as limited-angle, region-of-interest, sparse and exterior tomography?
Lecture 14 (March 17)
Choosing the regularization parameter. There are several methods, applicable to more or less of the different regularization approaches.
- L-curve method (mainly for Tikhonov)
- Morozov discrepancy principle (Tikhonov, has been extended to TV in recent research from Hong Kong)
- S-curve method (for sparsity-promoting regularization)
- Multi-resolution consistence method (for TV, see this file)
We will take a look of these methods and test some of them.
Lecture 13 (March 15)
Wavelet transform in dimensions 1 and 2. Large-scale sparsity-promoting reconstruction with Iterative Soft Thresholding Algorithm (ISTA).
Lecture 12 (Friday, March 3)
Planning the remaining lectures of the course. We used the "double team" brainstorming method and ended up with this set of suggestions:
The results of the planning are seen in the themes for lectures 13, 14 and 15 above.
Wednesday, March 1. The lecturer is sick (damn this flu season!).
Instead of the lecture, please read Section 3.2 of this note on the Haar wavelet transform: FourierSeries_Wavelets_v3.pdf
Study these Matlab routines and clarify to yourself how the codes and the text are connected.
Lecture 11 (February 24)
Recap of the construction of the tomographic "rectangle" phantom. Simulating tomographic data with no inverse crime, using higher-resolution simulation and interpolation. Applying
- Truncated SVD (TSVD),
- Tikhonov regularization, and
- Total Variation regularization
for the sparse-angle tomographic problem. Note that the Tikhonov regularization code above is written in a matrix-free fashion using the Conjugate Gradient (CG) method. We will talk about CG and other iterative solvers later in the course. The Total Variation computation is a two-dimensional adaptation of the quadratic programming formulation including inequality and equality constraints.
The files MyDS2.m and MyDScol.m are simple routines for downsampling an image by a factor of two.
See Sections 2.3, 4.4, 5.5, 6.2 and 5.2 in the book Mueller-Siltanen (2012).
THERE IS NO LECTURE on Wednesday, February 22. The lecturer is traveling.
Instead of the lecture, please run all the Matlab files given below related to Lecture 10 using this year's target signal AND with the photographic data collected in Lecture 9.
Lecture 10 (February 17)
Discussed were three forms of variational regularization: Tikhonov, generalized Tikhonov and Total Variation regularization. We focused on the one-dimensional convolution model and used a simple finite-difference matrix L in the regularization terms.
Matlab demonstrations were based on the routines developed in the 2015 Inverse Problems course. They are given below. I recommend you to run them through with this year's target and with the photographic data collected in Lecture 9.
Lecture 9 (February 15)
Photographic data with varying levels of blur (from defocusing of the lens) and noise (from increasing the ISO) as measurements.
We photographed printouts of these two pictures posted to the wall:
Further, we picked out rows from the edge image and analysed the pixels values in Matlab.
Lecture 8 (February 10)
Naive inversion for non-square matrices. More precisely, we studied the case when A is a k x n matrix with k>n and when all the singular values of A are strictly positive. In that case there is a unique least-squares solution f0 to the equation m = Af. Further, we have f0 = inv(A^T A) A^T m.
See Sections 4.1 and 5.2 in the book Mueller-Siltanen (2012).
Lecture 7 (February 8)
Building the computational model Af=m for the 2D tomography problem. That is,
-designing a phantom for test cases (we went for three rectangles in empty background in the unit square), called SqPhantom.m,
-constructing and saving the matrix A for a selection of sizes of f and collections of X-ray projection directions,
-calculating and saving the Singular Value Decomposition for the matrix A.
Also, we tried out Truncated Singular Value Decomposition (TSVD) for robust reconstruction.
Lecture 6 (February 3)
Recap of the trinity
(1) measurement data vector given by a device,
(2) continuum model of the physical measurement process, and
(3) computational model used in practical inversion.
The aim is to build two different example models, namely 1D deconvolution and 2D tomography. These two cases will be used throughout the course.
Construction of a matrix model m=Af using Matlab's radon.m routine. The matrix is constructed column by column by applying radon.m consecutively to unit vectors. While computationally inefficient, this approach provides an easy way for constructing system matrices for tomographic problems. The matrix A can then be studied by, for example, singular value decomposition.
Lecture 5 (February 1)
The most important is "Tomographic Imaging with Sparse Data."
Also, please read Section 2.3 of the textbook (Mueller&Siltanen 2012).
Lecture 4 (January 27)
Lecture 3 (January 25)
Quick look at the first noise-robust reconstruction method of the course: truncated singular value decomposition (TSVD).
We test TSVD with the simple convolution example from Lecture 2. Book chapters: 2.1, 4.2 and 4.1.1.
First observation about TSVD: the reconstruction is always a linear combination of singular vectors. Singular vectors are columns of the matrix V appearing in the singular value decomposition.
Second observation about TSVD: if we use only very few singular values, the reconstruction is not very accurate but extremely robust against noise. This means that the reconstruction does not change much form measurement to measurement even when there is large random noise added to the data.
Third observation about TSVD: it is not clear how to choose the truncation index r_\alpha in general. This is indeed a hard problem, and we will see some approaches to that later in the course.
Also, we started to discuss the difference between the continuum model for convolution (book chapter 2.1.1) and discrete model for convolution (book chapter 2.1.2).
Lecture 2 (January 20)
Discrete convolution between two vectors. How noise is amplified in naive inversion attempts. Book chapters: 2.1.
Lecture 1 (January 18)
Introduction to inverse problems. Practicalities about the course.
2 Materials and methods
Section 1 should briefly explain the topic in a way accessible to a non-expert.
Section 2 is for describing the data and the inversion methods used.
In section 3 the methods of Section 2 are applied to the data described in Section 2, and the results are reported with no interpretation; just facts and outcomes of computations are described.
Section 4 is the place for discussing the results and drawing conclusions.
The project is either about X-ray tomography or digital image processing. You can measure a dataset yourself in the Industrial Mathematics Laboratory.
X-ray topic: you can choose your own object to image in the X-ray lab. We can offer a range of tried-and-tested objects and tools to tailor them to your liking. Also, you can come up with your own idea of the measured object. The size of the object should not much exceed the size of an egg, and the chemical composition is important as X-ray attenuation contrast arises from electron densities in the materials. Please contact Alexander Meaney to find out if your object is good for imaging.
Photographic topic: take blurred or noisy photos of suitable targets.
In the project you are supposed to apply some of the regularization methods discussed in the course. Optimally, you should have an automatic method for choosing the regularization parameter and an automatic stopping criteria for the iteration. These are both difficult requirements, so have a simple approach as plan B if a more complicated approach does not work. Also, when choosing your objects of measurement, it's good to think about mathematical models of a priori information (piecewise constant, smooth, piecewise smooth) as it affects the choice of the regularizer.
First goal: the two first sections (Introduction and Materials and Methods) should be preliminary written in LaTeX, but not necessarily in poster format. The most important things to explain are:
- what kind of data to measure,
- what inversion method to apply for the reconstruction, and
- how to implement the computation.
The grade of the first goal represents 30% of the final grade of the project work. Please agree on a meeting time (in the period April 5-7) with the lecturer for reviewing and grading the first goal.
Second and final goal: poster session. The poster will be printed in size A1. You may create your own poster (from scratch), or you can use e.g. this template as a starting point and edit its layout, colors, fonts, etc. as much as you like.
The Matlab exercises will be checked individually at the exercise class.
Points for all completed exercises will be accumulated and part of the final grade.
Course feedback can be given at any point during the course. Click here.