README.md

Python tomography library

Introduction

This tiny Python library can be used to simulate and reconstruct X-ray tomography measurements. The library contains methods to calculate maximum a posteriori and conditional mean estimates out of sinograms of a given image.

Dependencies

• Python 3.6+
• Cython and C++ compiler with OPENMP support
• Numpy
• Scipy: Scipy sparse, Scipy signal, Scipy optimize
• PyWT
• Tqdm
• Skimage (scikit-image)
• Matplotlib
• H5py
• Pathlib
• Argparse

Usage

The main script can be used with or without command line arguments. When command line arguments are used, only one instance of the tomography class is created and only one reconstruction can be done. Alternatively one can type and call the desired calculations directly from the main.py, if multiple calculations are done and multiple instances of measurements are created.

Tomography class

The script's core is the tomography class, which is initialized by giving the target image (filename), size (targetsize), noise level (noise) and measurement angles (itheta). Inverse crimes can be avoided by simulating the sinogram by a different grid and reconstructing the image with a second one.
Parameter crimefree is False by default, dimbig and N_thetabig are the dimensions of the simulation sinogram. Preferrably they should be primes. The lhdev parameter refers to the sinogram measurement likelihood sigma. By default it's None, which means that the sigma is assumed to be proportional to the noise level: it's (sinogram_maxvalue * noise) ^ 2. If the noise level is 0, then the sigma is (sinogram_maxvalue * 0.01) ^ 2.

There are two ways to enter the measurement angles: user can enter either one integer (which then is the number of angles between 0 and 180 degress) or three integers, which will refer the first angle, last angle and the number of angles between them, respectively.

Example usage

For example, to initialize simulated sinogram measurement of the Shepp-Phantom image (64x64) with 59 angles and 5 percent noise without inverse crime, one shall run in the main.py:

theta = 50
t = tomography("shepp.png", 64, theta, 0.05, crimefree=True,commonprefix='/results/')

Then one can reconstruct the image from the sinogram with SCAM method plot the CM estimate and save the result details to a HDF5 file to the '/results/' subdirectory in the current directory.

ret = t.mwg_tv(2.5, 20000, 10000, thinning=10, mapstart=False, retim=False)
t.saveresult(ret)
plt.imshow(ret.result)
plt.show()

One can also obtain only a MAP estimate and just plot it:

ret2 = t.map_cauchy(0.01)
plt.imshow(ret2)
plt.show()

Compilation

The Cython files are compiled with the command

python3 setup.py build_ext --inplace

It will produce two compiled library files into the current directory, cyt and matrices. (i.e. in Linux cyt.so and matrices.so). Then the library is ready to use. Cyt-library contains functions for log-PDFs of different posteriors, their gradients and MCMC functions. Matrices-library contains functions for the Radon matrix operator and 2D Wavelet matrix operator construction.

Methods

After the class is initialized, the calculations itself can be run. The names of methods are rather self-descriptive. The methods of the tomography class beginning with with the word map refer to MAP estimates with different priors, mwg-starting methods refer to CM estimation by Metropolis-within-Gibbs (SCAM) and hmc-beginning methods refer to CM estimation by Hamiltonian Monte Carlo.

For all methods, the most important parameter is the prior's regularization parameter. Depending on the prior, it can be selected also so that the reconstructions remain approximately the same even if the resolution of the image is increased. The second important parameter is technical. With the retim parameter one can select what the methods actually return. The retim is by default True for all methods, which means that the methods return only the reconstructed image. If the parameter is False, an instance of container class is returned. The class can be feed as an argument to the saveresult method, which then saves the result to a hdf5 file for later analysis.

The container consists of attributes which are iterated through when saving to the result file. The class has only one method in addition to the init, finish. An instance of this class is supposed to be created when the calculation is started and the finish method is called when it's completed. This saves the execution time to one attribute.

However, everything interesting for most users is within the tomography class and its methods:

• tomography.__init__filename, targetsize=128, itheta=50, noise=0.0, commonprefix="", dimbig = 607, N_thetabig=421, crimefree=False,lhdev=None)

  • One might increase dimbig and N_thetabig from their default values, if the targetsize is near to 500.
  • commonprefix is a relative directory to the main.py, where the result files are saved (i.e. /results/ means $(PWD)/results/).

• map_tikhonov( alpha=1.0, order=1,maxiter=400,retim=True)

  • MAP function for Tikhonov regularization. Order means the order of the discrete derivative (1 or 2).
  • Maxiter parameter is the maximum number of minimization iteration rounds. Default value ensures that a resonable result is obtained within reasonable time. If the maxiter value is reached, one might increase it to ensure converge.

• map_tv( alpha=1.0, maxiter=400,retim=True)

  • MAP function for total variation regularization.

• map_cauchy(alpha=0.05, maxiter=400,retim=True)

  • MAP function for Cauchy difference prior. Default alpha value of 0.05 seems to have rather strong effect.

• map_wavelet(alpha=1.0, type='haar', maxiter=400,levels=None ,retim=True):

  • MAP function for total wavelet regularization. Default DWT level is None, which means that floor(log2(targetimagedim)-1) levels are used.

• hmcmc_tikhonov( alpha, M=100, Madapt=20, order=1,mapstart=False,thinning=1,retim=True)

  • HMC function for CM estimation with Tikhonov regularization. Since MAP and CM should converge to the same solution with Gaussian priors, this function is just for testing purposes and thus not interesting.
  • If retim is False, only the CM estimate result is returned, otherwise a container object is returned (which of course includes the CM, along with the chain).
  • Note that M and Madapt are almost certainly too small even for HMC, but enough to verify that the function works.
  • Mapstart decides if the MCMC chain is started from a MAP estimate rather than a random point.
  • Thinning is the reduction factor of MCMC samples, if the chain is wanted (retim=False). For HMC, thinning is usually not needed at all, since the chain is so short.

• hmcmc_tv(alpha, M=100, Madapt=20,mapstart=False,thinning=1,retim=True)

  • HMC function for CM estimation with TV regularization.

• hmcmc_cauchy(alpha, M=100, Madapt=20,thinning=1,mapstart=True,retim=True)

  • HMC function for CM estimation with Cauchy difference prior.

• hmcmc\wavelet( alpha, M=100, Madapt=20, type='haar',levels=None,mapstart=False,thinning=1,retim=True)

  • HMC function for CM estimation with Wavelet prior.

• mwg_tv(alpha, M=10000, Madapt=1000,mapstart=False,thinning=10,retim=True)

  • Metropolis-within-Gibbs (SCAM) function for CM estimation with Total variation prior.
  • With MwG, if one is willing to get the chain, thinning is usually required to save memory (from 10 to 500+, if a very long chain is wanted).

• mwg_cauchy( alpha, M=10000, Madapt=1000,mapstart=False,thinning=10,retim=True)

  • Metropolis-within-Gibbs (SCAM) function for CM estimation with Cauchy difference prior.

• mwg_wavelet( alpha, M=10000, Madapt=1000,type='haar',levels=None,mapstart=False,thinning=10,retim=True)

  • Metropolis-within-Gibbs (SCAM) function for CM estimation with Wavelet prior.

• sinogram()

  • Function, which plots the simulated sinogram without further actions.

• target()

  • Returns the target image of the simulation, which is resized to the targetsize.

• saveresult(result)

  • This is rather important function. It takes care of saving a given container result object to a HDF5 file.