Matrix-free X-ray tomography with sparse data

Matrix-free X-ray tomography with sparse data

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 solve simulated tomographic problems without constructing the measurement matrix A at all.
Furthermore, we collect data with sparse angular sampling, leading to extremely ill-posed reconstruction
problems.

Simulate tomographic data without inverse crime

This routine uses the radon.m function of Matlab's Image processing toolbox to create a sinogram:

XRsparseA_NoCrimeData_comp.m

As an example we simulated data at resolution 512x512 and adding a high level of noise 
(10%, achieved by setting noiselevel = 0.1; in the file XRsparseA_NoCrimeData_comp.m). 
We use 12 parallel-beam projection directions evenly distributed over 180 degrees.

The result looks like this (original Shepp-Logan phantom on the left, sinogram on the right): 

Compute reconstruction using filtered back-projection

Use this Matlab routine to compute filtered back-projection reconstructions:

XRsparseB_FBP_comp.m

As an example we simulated data at resolution 512x512 and adding a high level of noise 
(10%, achieved by setting noiselevel = 0.1; in the file XRsparseA_NoCrimeData_comp.m). 
Running the file XRsparseB_FBP_comp.m gives the following output (original Shepp-Logan
phantom at left, FBP reconstruction on the right): 

We remark that the filtered back-projection (FBP) methodology is not at all designed for sparsely
sampled data considered here. FBP works very well with denser angular sampling, whenever that
is available. 

Also, probably FBP could be optimized to give better results even with the sparse data used here.
However, that is outside the scope of this treatment. 

Compute reconstruction using matrix-free iterative Tikhonov regularization

XRsparseC_Tikhonov_comp.m, XRsparseC_Tikhonov_plot.m

As an example we simulated data at resolution 512x512 and adding a high level of noise 
(10%, achieved by setting noiselevel = 0.1; in the file XRsparseA_NoCrimeData_comp.m). 
Then we set the following parameter values in the file called XRsparseC_Tikhonov_comp.m: 

K = 600;         
alpha = 300;

The result looks like this (original Shepp-Logan phantom at left, Tikhonov reconstruction on the right):

Compute reconstruction using matrix-free iterative (approximate) total variation regularization

Here are the main reconstruction and plot routines:

XRsparseD_aTV_comp.m, XRsparseD_aTV_plot.m

You will also need these files:

XR_aTV_feval.m, XR_aTV_fgrad.m, XR_aTV_grad.m, XR_aTV.m, XR_misfit_grad.m, XR_misfit.m

As an example we simulated data at resolution 512x512 and adding a high level of noise
(10%, achieved by setting noiselevel = 0.1; in the file XRsparseA_NoCrimeData_comp.m).
Then we set the following parameter values in the file called XRsparseD_aTV_comp.m: 

alpha = 30000;
MAXITER = 12000;               
beta    = .000001;

The result looks like this (original Shepp-Logan phantom at left, approximate total variation reconstruction on the right):