X-ray tomography with matrices

X-ray tomography with matrices

This page contains the computational Matlab files related to the book
 Linear and Nonlinear Inverse Problems with Practical Applications
 written by Jennifer Mueller and Samuli Siltanen and published by SIAM in 2012.

You can order the book at the SIAM webshop.

Go to master page

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:

XRMA_matrix_comp.m

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:

XRMF_truncSVD_comp.m

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