The SIRT Algorithm

This article follows on “Tomography, Part 4: Algebra!” of the series on tomography. Reread the complete series starting from part 1 if you need a refresher. As a practical example of an algebraic technique, I present the SIRT (Simultaneous Iterative Reconstruction Technique) algorithm. As in the last article of the series on tomography, I start from the system of linear equations,

\[{\bf Ax}={\bf b},\]

where \({\bf x}\) represents the image, \({\bf b}\) represents the projections, and \({\bf A}\) represents the scanning process (see the mentioned article for details).

A system like the one shown is typically solved by minimizing some norm \(\|{\bf Ax}-{\bf b}\|\). The SIRT algorithm is one of many methods for doing just that.

Definition

First recall that \({\bf A}\) represents the action of the scanner, also known as the forward projection. Each row of \({\bf A}\) contains the coefficients of an equation that corresponds to a single ray. It describes how the pixels are combined into ray sums. The transposed matrix, \({\bf A}^\mathrm{T}\) back projects the projection images onto the reconstruction area. Given a ray sum, it describes which pixels are hit by that ray.

SIRT alternates forward and back projections. Its update equation is

\[{\bf x}^{(t+1)}={\bf x}^{(t)}+{\bf CA}^\mathrm{T}{\bf R}({\bf b}-{\bf Ax}^{(t)}),\]

where \({\bf C}\) and \({\bf R}\) are diagonal matrices that contain the inverse of the sum of the columns and rows of the system matrix, so \(c_{jj}=1/\sum_i{a_{ij}}\) and \(r_{ii}=1/\sum_j{a_{ij}}\). These matrices compensate for the number of rays that hit each pixel and the number of pixels that are hit by each ray.

Iterations start with \({\bf x}^{(0)}={\bf 0}\). The update equation then does the following things.

1. The current reconstruction \({\bf x}^{(t)}\) is forward projected: \({\bf Ax}^{(t)}\).
2. The result is subtracted from the original projections: \({\bf b}-{\bf Ax}^{(t)}\).
3. This difference is then back projected. In essence this is done by multiplying with \({\bf A}^\mathrm{T}\), but weighted with \({\bf C}\) and \({\bf R}\).
4. This results in the correction factor \({\bf CA}^\mathrm{T}{\bf R}({\bf b}-{\bf Ax}^{(t)})\).
5. The correction factor is then added to the current reconstruction, and the whole process is repeated from step 1.

And that’s it. It’s not more complicated than that.

Implementation

For demonstration purposes, you can implement the update equation directly in MATLAB or Python. Below is an implementation in MATLAB, in which I have assumed that the system matrix \({\bf A}\) is known.

% Input: sparse system matrix A, data b.
% Output: SIRT reconstruction x.
x = zeros(d * d, 1);
[rows cols] = size(A);
C = sparse(1 : cols, 1 : cols, 1 ./ sum(A));
R = sparse(1 : rows, 1 : rows, 1 ./ sum(A'));
CATR = C * A' * R;
for i = 1 : 100
  x = x + CATR * (b - A * x);
end

This code is only useful for very small experiments, but it has the big advantage of being very close to the mathematical expression. For use with practical datasets, you can use the ASTRA Tomography Toolbox, which contains a GPU-accelerated implementation. That Toolbox can also provide you with the system matrix of any geometry that you define in it.