Summary: Tomography means reconstructing the internal structure of a physical body using X-ray images of the body taken from different directions. Mathematically, the problem is to recover a non-negative function f(x) from a collection of line integrals of f. Here we show how to solve the tomography problem using Total Variation (TV) regularization on experimental data. TV is an edge-preserving reconstruction method for ill-posed inverse problems that favors sharp piecewise constant reconstructions. TV was first introduced by Rudin, Osher and Fatemi in 1992 [1].
Authors: Kristian Bredies, Tatiana Bubba and Samuli Siltanen.
Software: Tested on Matlab 2018a.
Data: The open data set lotusroot.
Download package as zip file: CodesTV
Literature: L. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60, 259–268 (1992).
Please read also the post Total variation regularization for X-ray tomography for an introduction on TV tested on simulated data. A summary of the experimental data can be found on arXiv.
The usual Total Variation (TV) seminorm [1] is defined as
TV is usually a suitable model for images with edges and have become a very popular tool in image processing and inverse problems due to peculiar features that cannot be realized with smooth regularizations [2].
In the following we compare the reconstructions obtained by minimizing the penalty functional:for different values of the regularization parameter α. Here, A is the measurement matrix, m the experimental data, f a discretization of the object to reconstruct, and ι≥0 the indicator function of the nonnegative orthant. The discretization of TV formulas is carried out by finite differences as described in details in [3].
For the solution of the minimization problem above, any kind of numerical algorithm for the solution of convex-concave saddle point problems can be used. The algorithm implemented in the dowloadable package above is a primal-dual ascent-descent method with primal extragradient, as described in the famous 2011 paper by Chambolle and Pock [4]. The main bulk of the code is in tomo_tv.m. The script recon_tv.m can be used to set up the experiment: loading the data and setting up the values for different parameters (the regularization parameter and the maximum number of iterations).
In the following computational examples we use the 256×256 lotus root dataset obtainable from here (documentation in arXiv). This dataset uses 120 projection directions. We set the maximum number of iterations equal to 10000 and the regularization parameter α = 10−5 (see Figure 1).
Figure 1: Reconstructions with α = 10−5 .
As always, the value for the regularization parameter is chosen to balance the data mismatch term and the prior information given by the TV penalty. We can play around with the value of α (still keeping the maximum number of iterations equal to 10000) to see how this affects the reconstructions (see Figures 2).
Figure 2: Reconstructions with α = 10−4 (top) and α = 10−6 (bottom).
Unfortunately, sometimes it can also happen that TV reconstructs undesired edges: this artifact is called staircasing effect. This is due to the fact that the model assumption for TV is that an image is piecewise constant up to a discontinuity set. However, natural images are often piecewise smooth due to shading, for instance. Stay tuned for the next post to see how this can be overcome.
References
[1] L. Rudin, S. Osher and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60, 259–268 (1992).
[2] M. Burger and S. Osher, A guide to the TV zoo (chapter 1 in M. Burger and S. Osher, Level-Set and PDE-based Reconstruction Methods, Springer, 2013).
[3] K. Bredies, Recovering piecewise smooth multichannel images by minimization of convex functionals with Total Generalized Variation penalty, Lecture Notes in Computer Science 8293, 44–77 (2014).
[4] A. Chambolle and T. Pock, A first-order primal-dual algorithm for convex problems with applications to imaging, Journal of Mathematical Imaging and Vision 40, 120–145 (2011).
