Pseudo-iterative method : SIRT-Filter reconstruction¶
Dependencies¶
For the filters computation, this feature needs the ASTRA toolbox and its Python wrapper to be installed. If this is a problem for you, the filters can be pre-computed outside PyHST with a Python script on another machine. On the ESRF cluster the only small problem is that we have not yet a debian package for astra (which is used to calculate the sirt filter )
Therefore if you want to test this new feature with pyhst2_unstable you have to add these two paths to your .bashrc
for DEBIAN7
export PYTHONPATH=/scisoft/ESRF_sw/astra/debian7/python/
LD_LIBRARY_PATH=/scisoft/ESRF_sw/astra/debian7/lib:$LD_LIBRARY_PATH
and for DEBIAN8
export PYTHONPATH=/scisoft/ESRF_sw/astra/debian8/python/
LD_LIBRARY_PATH=/scisoft/ESRF_sw/astra/debian8/lib:$LD_LIBRARY_PATH
Principle¶
In parallel geometry, the SIRT iterative algorithm can be approximated by a one-step filtering operation (see Reference).
Assume you want to run \(N\) iterations of the SIRT algorithm on a volume. First, a filter is computed. This filter only depends on the geometry of the dataset, namely :
The slice dimensions \((n_x, n_y)\)
The horizontal dimension of the detector \(n_d\)
The number of projections \(n_A\)
The rotation axis position
This filter is computed once for one slice, and is applied to the all slices of the volume. Moreover, if another volume has the same geometry, the same filter can be used. PyHST has a saving mechanism that enables to use a pre-computed filter (you can therefore have a database of available filters). In practice, this means that the SIRT algorithm can be replaced by another equivalent method, with the cost of a simple Filtered-Backprojection.
The filters are saved in the form sirtfilter_nx_ny_nA_niter.h5, as a HDF5 file. More information is available in the attributes of the HDF5 file.
Note : if the Python module h5Py is not available, the filters are stored as a compressed numpy file (.npz).
Syntax¶
To use this reconstruction type, the following parameters are relevant :
ITERATIVE_CORRECTIONS=0 # Disable all other methods
DENOISING_TYPE = 9 # Parameter enabling SIRT-Filter
SF_ITERATIONS = 300 # Number of SIRT iterations
SF_FILTER_FILE = /path/to/the/filters/ # Path where the filters are stored/loaded
SF_LAMBDA = 0.0 # Tikhonov regularization
Note that it is important to set ITERATIVE_CORRECTIONS=0. Also, the parameter SF_FILTER_FILE is actually a directory, not a file.
Tikhonov regularization¶
It is well known that the SIRT algorithm does not provide stable solutions - that is, if the number of iterations if too high, the noise level in the resulting slice is likely to be largely amplified.
A regularization term can be added to stabilize the solution.
An example is the Tikhonov regularization \(\lambda \norm{x}_2^2\).
The parameter SF_LAMBDA corresponds to \(\lambda\), the regularization weigth.
Reference¶
This method is based on this work :
Pelt, D. M., & Batenburg, K. J. (2015). Accurately approximating algebraic tomographic reconstruction by filtered backprojection. Proceedings of the 2015 International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine.
Mathematical background¶
SIRT is an algorithm solving the least-squares problem
Where \(d\) is the acquired data, and \(x\) is the slice. The iteration \(k\) of a simple gradient descent algorithm is :
which means that the iteration \(k+1\) can be expressed linearly as a function of the iteration \(k\) :
where \(A = \left(\id - \gamma P^T P\right)\) and \(b = \gamma P^T d\). The slice at iteration \(n\) can then be directly obtained :
The basic idea is to pre-compute the operator \(\mathcal{A} = \left(\displaystyle\sum_{k=0}^{n-1} A^k \right)\). If this can be done, then the whole SIRT algorithm could be performed in one step :
which is a “backproject-then-filter” method. The operator \(\mathcal{A}\) only depends on the geometry.
In general, if \(N^2\) is the size of a slice, then \(\mathcal{A}\) would have a size \(N^2 \times N^2 = N^4\), which cannot be stored (275 GB would be necessary for \(512 \times 512\) slice, and 70 TB for a \(2048 \times 2048\) slice). However, in parallel geometry, the operator \(P^T P\) involved in \(\mathcal{A}\) has a very special structure : it is nearly a Toeplitz matrix, meaning that this operator can be (nearly) represented by a convolution kernel of size \(\sqrt{N^2 \times N^2} = N^2\), the size of a simple image.
The filter kernel is computed from a Dirac delta function. Then, the kernel is projected in the sinogram domain, and its Fourier transform is computed. The real part of the FT is stored and will be used as the filter. To reconstruct a given sinogram, the FT of the sinogram is computed, then multiplied with the filter. After computing the IFT, the resulting sinogram is back-projected. The result is the reconstructed slice.