 |
Linear and Nonlinear Inverse Problems with Practical ApplicationsJennifer L. Mueller and Samuli Siltanen |
X-ray Tomography with Matrices
This page contains the computational Matlab files related to the book, which you can purchase in the SIAM Bookstore.
Here we construct the measurement matrix A related to tomographic imaging and solve the inverse problem using several regularized methods. These routines make heavy use of radon.m and phantom.m files that are available only in Matlab's Image processing toolbox, not in the basic version of Matlab.
Constructing the measurement matrix A
This Matlab routine builds the tomographic measurement matrix for NxN images and N uniformly distributed parallel-beam projections:
You can choose your favorite value of N inside the above m-file. Note that taking N greater than 64 will choke most personal computers either here or later when computing the singular value decomposition of A! The examples on this page are meant to be studied at low dimensions, such as N between 16 and 50, say.
The matrix is saved to your working directory as a .mat file. The number N is part of the filename, as you can see in these example files:
RadonMatrix16.mat
RadonMatrix32.mat
Creating data that contains inverse crime
You can try out naive inversion involving inverse crime with the following routines:
XRMB_naive_comp.m
XRMB_naive_plot.m
Creating data without inverse crime
These routines interpolate the data from a higher-dimensional model.
XRMC_NoCrimeData_comp.m
XRMC_NoCrimeData_plot.m
You can go ahead and see how naive inversion fails again:
XRMD_naive_comp.m
XRMD_naive_plot.m
Regularized inversion using truncated SVD
Compute the singular value decomposition:
XRME_SVD_comp.m
XRME_SVD_plot.m
Now we can use the truncated SVD to achieve robust reconstructions:
Tikhonov regularization
Compute Tikhonov regularized reconstructions with these routines:
XRMG_Tikhonov_comp.m
XRMG_Tikhonov_plot.m
You can experiment by varying the regularization parameter.
Approximate total variation regularization
Here we reconstruct the attenuation coefficient using total variation regularization. We smooth the non-differentiable corner of the absolute value function appearing in the total variation norm, resulting in a smooth objective functional to be minimized. The optimization method we use is the Barzilai-Borwein method.
XRMH_aTV_comp.m
XRMH_aTV_plot.m.
You will need all of the files below to run the above reconstruction routine:
XRMH_aTV_feval.m
XRMH_aTV_fgrad.m
XRMH_aTV_grad.m
XRMH_aTV.m
XRMH_misfit_grad.m
XRMH_misfit.m
