Inversio-ongelmat, kevät 2012
Inversio-ongelmat, kevät 2012
Inverse problems, spring 2012
Finnish summary
Inversio-ongelmien tutkimus on sekä puhtaan että sovelletun matematiikan aktiivinen
osa-alue, joka on Matematiikan ja tilastotieteen laitoksella edustettuna kolmen professorin,
lukuisten tutkijoiden ja jatko-opiskelijoiden sekä Suomen Akatemian huippuyksikön voimin.
Alan tutkimuksessa tavoitteena on saada tuntemattomasta kohteesta tietoa epäsuorien mittausten avulla.
Inversio-ongelmilla on sovelluksia monilla eri tieteen ja tekniikan aloilla.
Näitä ovat mm. signaalinkäsittely, rakenteita tuhoamaton testaus, fysiikka ja astrofysiikka,
muodon optimointi, lääketieteellinen kuvantaminen ja potilaan elintoimintojen seuraaminen,
solubiologia, geofysikaalinen kuvantaminen sekä optiohinnottelu osakemarkkinoilla.
Esimerkkejä inversio-ongelmista:
-Epätarkan valokuvan terävöittäminen
-Kolmiulotteinen röntgentomografia
-Maan rakenteiden kuvantaminen maanjäristysaaltojen kulkuaikojen avulla
-Halkeamien etsiminen rakenteista
-Öljyn ja mineraalien etsintä
-Saasteiden maanalaisen leviämisen havainnointi
-Asteroidien muodon selvittäminen kirkkauden muutosten avulla
Yhteistä kaikille näille ongelmille on niiden herkkyys
mittausvirheille sekä se, että tutkittavasta kohteesta saadaan
vain epäsuoria ja kohinaisia mittauksia.
Keväällä järjestetään laaja inversio-ongelmien kurssi, jonka
voi sisällyttää sekä yleisen että sovelletun matematiikan syventäviin opintoihin.
Kurssin luennot koostuvat kahdesta osasta:
-Analyyttiset inversio-ongelmat (prof. Matti Lassas)
-Laskennalliset inversio-ongelmat (prof. Samuli Siltanen).
Luentojen laajuus on 5+5=10 op.
Luentojen lisäksi kurssi sisältää aiempaa "Matematiikan sovellusprojektit"-kurssia
vastaavan projektityön. Tämä tehdään ohjattuna parityönä ja sen laajuus on 5 op.
Yhteensä kurssin laajuus on siis 15 op.
English summary
Inverse problems research is an active area of both pure and applied mathematics.
At the Department of Mathematics and Statistics the field is represented by three professors,
several researchers and graduate students and a Centre of Excellence of Academy of Finland.
Inverse problems are about interpreting indirect measurements, and they have
applications to many areas of science and technology. Examples include signal processing,
non-destructive testing, astrophysics, medical imaging and geophysical prospecting.
The lectures of the course consist of two parts:
-Analytic inverse problems (Prof. Matti Lassas)
-Computational inverse problems (Prof. Samuli Siltanen).
Together the lecture parts make up 5+5=10 credit units.
In addition to lectures the course involves a project work. It is done in teams of two and
gives 5 credit units to each student.
The course is in total 15 credit units.
Lecturers
Matti Lassas and Samuli Siltanen
Credit units
15 op.
Type
Advanced studies (syventävä opintojakso)
Exam
Wednesday 16.5., 10-12, B120.
Prerequisites
Recommended: basics of Fourier transform, measure theory, integration and linear algebra
(for example courses Fourier-muunnoksen perusteet, Mitta- ja integrointiteoria, Lineaarialgebra 1&2).
Experience with Matlab is helpful in the computational part of the course, but not necessary.
Tutorials for learning basics of Matlab are available for example at the following locations:
Lectures
First lecture will be held Wednesday 18.1.2012
Period III: Lectures and exercises, weekly schedule can be found below.
Tuesday 10-12, Exactum hall B120,
Wednesday 10-12, Exactum hall B120,
Friday 12-14, Exactum hall D123.
Wed 18.1. Matti Lassas: Introduction and applications of inverse problems
Fri 20.1. Matti Lassas: General questions studied in inverse problems. Definition of X-ray and Radon-transforms
Tue 24.1. Samuli Siltanen: One-dimensional convolution, both continuous and discrete. Book: chapter 1 and sections 2.1.1 and 2.1.2.
Wed 25.1. Samuli Siltanen: One-dimensional deconvolution: naive inversion, ill-posedness and singular value decomposition. Book sections 2.1.3 and 2.1.4 and Chapter 3.
Fri 27.1. Matti Lassas: Fourier-series.
Tue 31.1. Matti Lassas Fourier-transform (basics)
Wed 1.2. Samuli Siltanen: Truncated singular value decomposition as a regularization method. Book Chapters 3 and 4.
Fri 3.2. Samuli Siltanen: X-ray tomography. Book Chapters 2-4. Slides:
Tue 7.2. Matti Lassas: Fourier-transform (inverse transform, mapping properties in Schwartz class and L2 spaces)
Wed 8.2. Matti Lassas Fourier-slicing formula for Radon transform.
Fri 10.2. Samuli Siltanen: Tikhonov regularization. Book Chapter 5.
Tue 14.2. Samuli Siltanen: Tikhonov regularization. Book Chapter 5.
Wed 15.2. Matti Lassas Properties of temperated distributions and Fourier transform for distributions
Tue 28.2. Matti Lassas Sobolev spaces on a Mobius band and properties of Radon transform
Wed 29.2.-Fri 9.3: No Lectures
Tue 13.3. Matti Lassas: Filtered backprojection method
Wed 14.3. Matti Lassas: Singularities in X-ray imaging and inverse problems for the wave equation
Fri 16.3. Matti Lassas: Inverse problems for the wave equation
Later part of period IV: Project work.
Materials
Material for the analytical part of the course will be uploaded here. The material
contains:
Exercises
1.Exercises of the analytical part
Deadline Fri 17.2.
1. Matlab-exercise, deadline Fri 10.2. Please return at the lecture the printouts as detailed below.
Start by running the following files below related to one-dimensional deconvolution:
DC1_cont_data_comp.m, DC2_discretedata_comp.m, DC3_naive_plot.m, DC4_truncSVD_comp.m.
Next you are going to modify those files.
First make the parameter "a" larger inside the file DC1_cont_data_comp.m. This makes the blurring effect
of the convolution much stronger, and the inversion task becomes harder. Does the naive inversion with
inverse crime fail with large enough "a"? If so, find the smallest "a" for which this happens and print out
a picture of the data and of the failed reconstruction. Does the discretization parameter "n" in file
DC2_discretedata_comp.m have some effect on the above? How do the singular values behave in
file DC4_truncSVD_comp.m for a choice of "a" for which the inverse-crime reconstruction fails?
Next make the parameter "a" smaller in file DC1_cont_data_comp.m, and choose a quite fine discretization
in file DC2_discretedata_comp.m, for example take n=256. How do the singular values behave when "a"
is made smaller? Print out three plots of singular values with decreasing "a" to show what happens.
Extra challenge for more advanced students (not mandatory): replace the trapezoidal rule on line 42
in the file DC1_cont_data_comp.m by Simpson's rule. Compare the accuracy of the result.
If you are familiar with Matlab, you can do the exercise by yourself. If you need help, please contact
Esa.Niemi(at)helsinki.fi; he will organize a tutorial session. The tutorial session takes place on Tuesday 7th Feb, 16-18, hall C128.
2. Theoretical exercise, deadline Wed 1.2.
Take f
to be the function having value 1 in the interval 1/4 < x < 3/4 and zero otherwise.
Take a=0.1 and let the PSF to be the function having value 1 in the interval -a < x < a.
Compute the convolution between the PSF and function f by hand directly from Definition 2.1.1 of the book.
3. Matlab exercise session, Tue 21.2. 10-12, hall C128
Matlab routines: , , ,
Some tomography matrices computed with (for those who do not have radon.m available).
, , , .
Also, here are some results created by the routine :
, , , .
The idea of this exercise is to illustrate some basic functions and operations in Matlab by implementing the TSVD regularized solution for the linear equation Ax=mncn (given the matrix A and noisy measurement mncn). Also, one should note that matrices, and thus SVDs, are not practical in 2D (or 3D) problems; computing the SVD is computationally too expensive.
Download the routines XRM1_matrix_comp.m, XRM3_NoCrimeData_comp.m and XRM6_truncSVD_comp.m. Ensure that the resolution is N=32 in the first two of the files.
Run first the routine XRM1_matrix_comp (this may take some time) and then XRM3_NoCrimeData_comp to create the matrix A and a noisy measurement mncn for the model Ax=mncn of x-ray tomography.
Now you should compute the TSVD solution of Ax=mncn by using M singular values (vectors). To do this, add necessary code after line 25 to the routine XRM6_truncSVD_comp.m; the code should compute the TSVD solution and save it to variable x. [Hint: use Matlab functions svd.m, size.m and basic matrix/vector operations.]
Having implemented the TSVD code to XRM6_truncSVD_comp.m
a) compute the TSVD solutions using 1, 10, 100, 500 and 700 singular vectors
b) what happens to the reconstruction with 700 singular vectors if the noiselevel (in XRM3_NoCrimeData_comp.m) is increased to 5% or to 10%?
(ii) In this exercise we compute the regularized solution by using Tikhonov regularization and stacked form. Moreover, we apply Morozov's discrepancy principle for choosing the regularization parameter. Notice that we still use matrices -> not very efficient computationally.
Download the routine XRM9_TikhStacked_comp.m (in addition to XRM1_matrix_comp.m, XRM3_NoCrimeData_comp.m). Ensure that the resolution is N=32 in files XRM1_matrix_comp.m and XRM3_NoCrimeData_comp.m. Set the noiselevel to 1% in XRM3_NoCrimeData_comp.m.
Run first the routine XRM1_matrix_comp (this may take some time) and then XRM3_NoCrimeData_comp to create the matrix A and a noisy measurement mncn for the model Ax=mncn of x-ray tomography.
Now you should compute the regularized solution of Ax=mncn by using Tikhonov regularization and stacked form. To do this, add necessary code after line 25 to the routine XRM9_TikhStacked_comp.m; the code should form the matrix and the right-side vector for the "stacked form" equation and solve it by the standard Matlab backslash operation.
Questions:
a) what size(s) of regularization parameter alpha give best reconstructions (by visual inspection)?
b) apply Morozov's discrepancy principle for choosing the regularization parameter alpha. What size of alpha does the principle suggest?
4. Matlab exercise session, Wed 22.2. 10-12, hall C128
Matlab routines: , .
Let us study the 1D heat conduction problem. First we consider the forward problem. In particular, we will see how the differences and errors in the initial temperature distributions get ``smoothed out'' as time proceeds.
Run the file HeatCondSimulator.m. Note that it plots two temperature distributions at some instant with ``pause'' mode: to see the distribution at the next time step, press some key. Repeat the same test but change the piecewise linear initial temperature distribution to a highly noisy version of the distribution 10*sin(2x) by uncommenting line 46 in HeatCondSimulator.m. What do you observe? Does the results tell something about the possible ill-posedness (instability) of the corresponding inverse problem?
(ii) Next let us try to solve the inverse heat conductivity problem; i.e, assume we are given the temperature distribution at time Tmax>0 and we should determine the initial temperature distribution. We do this for two different values of Tmax: Tmax=0.1 and Tmax=0.2. (Which one of the cases is expected to be more difficult?)
Ensure that Tmax = 0.1 and noiselevel = 0.0 in file InverseHeatCondNaive_withCrime.m and run it. Is the reconstruction good? Is there a reason to suspect an inverse crime?
Increase next the noiselevel to 0.0001 and compute the reconstruction again. What do you observe in terms of the inversion's sensitivity to noise? Increase now Tmax to 0.2 and set noiselevel back to zero. Does the naive inversion still succeed with the aid of an inverse crime?
5. Theoretical exercise, deadline Fri 24.2.
Matlab resources
One-dimensional convolution and deconvolution
These files are used by the other routines: constructs convolution matrices, defines the point spread function and plots it, and finally defines our basic signal we want to recover.
The first step is to study continuous convolution using highly accurate numerical integration. The file computes the approximation (you can modify the width of the PSF by changing parameter a), and the file plots the result.
After running the above "DC1"-files, you can simulate measurement data using . The file constructs three types of data: data with inverse crime for demonstrational purposes, very low-noise data without inverse crime, and data without inverse crime contaminated by white noise. You can change the noise amplitude by changing the parameter sigma. The data can be plotted using .
Naive inversion (applying the inverse of A to data) is demonstrated with file . This approach fails with ill-posed inverse problems.
The first regularized inversion method we study is truncated singular value decomposition, demonstrated by .
The second regularized inversion method is Tikhonov regularization, demonstrated by and, in a generalized form, by .
Here are Matlab resources for computing total variation regularized reconstructions using quadratic programming: . Note that Optimization toolbox is needed as the routine calls the quadprog.m function. Furthermore, you can study sparsity-based choice of regularization parameter with the files and .
Project work
The idea is to study an inverse problem both theoretically and
computationally in teams of two students. The end product is a
written report consisting of 10-20 pages. We recommend the classical table of contents for structuring the document:
1 Introduction
2 Materials and methods
3 Results
4 Discussion
Section 2 is for describing the data and the inversion methods used. In section 3 those methods are applied to the data 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.
First goal: two first sections should be preliminary written. Student A describes the simulation or collection of data, and Student B describes the reconstruction method that will be used. At this point no numerical computations are not necessary. The draft of project work presented at the first deadline will be graded, and the grade represents 30% of the final grade of the project work.
Second and final goal: report is complete. The idea is that Student B simulates or collects the data, and Student B implements the reconstruction.
Topics for project work:
1. Stability and continuity of Radon transform.
The work concentrates on understanding correct function spaces to study Radon transform and its inverse and possibly considering numerical examples. References: Bal's lecture notes, Natterer's book.
2. Resistor networks.
Study of electrical currents in networks of electrical resistors. The idea is to read some related literature (Curtis) and see how the problems can be written in matrix form. This is a topic in discrete mathematics related to medical imaging (electrical impedance tomography).
3. Invisibility.
4. Inverse problems and Riemannian geometry.
Understanding how geodesic and travel times are related or understanding geometric formulation of inverse conductivity problem. The work is related to invisibility cloaking.
5. Inverse spectral problem or non-linear inverse problem for wave equation.
The work is related to a classical question "Can you hear the shape of a drum?". In the project work one studies one-dimensional spectral problem or one dimensional wave equation and the goal is to understand some of the basic methods how the wave speed can be reconstructed using boundary measurements.
References: 1: Inverse boundary spectral problems (Katchalov-Kurylev-Lassas (Chapter 1), 2: An introduction to inverse scattering and inverse spectral problems (Chadan et al), Chadan's chapter, and 3: An introduction to the mathematical theory of inverse problems (Kirsch, Andreas, Applied Mathematical Sciences, 120. Springer-Verlag, New York, 1996), Chapter 4.
References 1-2 can be found from Kumpulan campus library and reference is provided by the lecturers. More material can be found from American Mathematical Societys web page Mathscinet by searching with words "inverse spectral". In the work one aims to write a report which deals either on the Boundary Control method (Reference 1) or the Gelfand-Levitan method (References 2 and 3).
6. Deblurring a misfocused photograph.
Take a printed paper with some text and an image; a newspaper page will do. Make sure that the page contains a small isolated black dot. It's a good idea to take a sharp picture for comparison and a few misfocused versions with different amount of blur. Also, take photographs with different sensitivities (ISO settings in the camera), resulting in various amounts of measurement noise. Read the photograph into Matlab using the "imread" command and pick out only one of the three color components (red, green or blue). Then you have a grayscale pixel image that can be viewed as a discrete approximation of a real-valued light intensity distribution on the paper. You can read off the shape of the point spread function from the image of the isolated black dot. Apply Tikhonov regularization or (approximate) total variation regularization and find out if you can make the blurred text readable. Read the photograph into Matlab using the "imread" command and pick out only one of the three color components (red, green or blue). Then you have a grayscale pixel image that can be viewed as a discrete approximation of a real-valued light intensity distribution on the paper. You can read off the shape of the point spread function from the image of the isolated black dot. Apply two regularization methods:
-Tikhonov regularization with the large-scale conjugate gradient approach,
-(approximate) L1 norm regularization by modifying the total variation code available,
and find out if you can make the blurred text readable.
7. Statistical inversion and Monte Carlo Markov Chain methods
Statistical (Bayesian) inversion provides efficient means for incorporating a priori information into reconstructions. In this project work, one-dimensional deconvolution data is inverted using smoothness requirement (as in generalized Tikhonov regularization with derivative penalty) together with enforcing the information that the signal takes values between zero and one.
8. Wavelets and Besov space regularization
Implement the Besov space regularization method with parameters s=p=q=1 for the one-dimensional deconvolution problem using approximate absolute value function and Barzilai-Borwein minimization. This can be tested on the 1D deconvolution example with a modest size of n, for example n=32. The Haar wavelet transform can be used, and it can simply be written as a matrix.
9. Total variation regularization for tomography
In the lecture material, total variation regularization was implemented for the one-dimensional deconvolution problem using the quadratic programming approach. The aim of this project work is to implement the same approach for the two-dimensional X-ray tomography application.
Registration
Forgot to register? What to do?.