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.

- The current reconstruction \({\bf x}^{(t)}\) is forward projected: \({\bf Ax}^{(t)}\).
- The result is subtracted from the original projections: \({\bf b}-{\bf Ax}^{(t)}\).
- 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}\).
- This results in the correction factor \({\bf CA}^\mathrm{T}{\bf R}({\bf b}-{\bf Ax}^{(t)})\).
- 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.

Could you please comment on the C and R matrix being related to the Gramian matrices. Is this so? Many thanks !

I don't think that there is any connection between the \({\bf C}\) and \({\bf R}\) matrices and the Gramian matrix… Do you have any indication that there should be one?

## Add new comment