Tomography, Part 4: Algebra!

This is the last of four articles on tomography (also read the first one, the second one, and the third one).

In the first three articles, I’ve focused on analytic techniques, in which the reconstruction problem is attacked using mathematical analysis. When the step to real world scanners is made, the problem is discretized. Both the projections and the reconstruction area are divided into pixels, and a numerical approximation of the true mathematical technique is introduced.

The algebraic techniques assume a discretized problem from the very beginning, and solve the reconstruction problem in a completely different way. In the algebraic approach, the reconstruction problem is represented by a system of linear equations. The variables (unknowns) are the pixels of the reconstruction area, and the right hand side of the equations is the projection data. A system of linear equations with \(n\) unknowns and \(m\) equations looks like

\[\begin{alignat}{7}
a_{11} x_1 &&\; + \;&& a_{12} x_2 &&\; + \cdots + \;&& a_{1n} x_n &&\; = \;&&& b_1 \\
a_{21} x_1 &&\; + \;&& a_{22} x_2 &&\; + \cdots + \;&& a_{2n} x_n &&\; = \;&&& b_2 \\
\vdots\;\;\; && && \vdots\;\;\; && && \vdots\;\;\; && &&& \;\vdots \\
a_{m1} x_1 &&\; + \;&& a_{m2} x_2 &&\; + \cdots + \;&& a_{mn} x_n &&\; = \;&&& b_m, \\
\end{alignat}\]

where the \(x_j\)’s (\(j=1,\ldots,n\)) represent the image, and the \(b_i\)’s (\(i=1,\ldots,m\)) represent the projection data. Using matrix notation, the linear system can be written more succinctly as

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

where \({\bf x}={(x_1\cdots x_n)}^\mathrm{T}\) and \({\bf b}={(b_1\cdots b_m)}^\mathrm{T}\). \(\mathrm{T}\) is the transpose operator, since we need \({\bf x}\) and \({\bf b}\) to be column vectors. So, you need to list all pixels of the reconstruction area in a single column vector \({\bf x}\), and all data points of the projections in a single column vector \({\bf b}\).

The System Matrix

What should the \(m\times n\) matrix \({\bf A}\) look like? This system matrix specifies exactly what the relation between the scanned object and the projections is. In other words, it is the mathematical representation of the scanner with which the projections were taken. Each row of the matrix (each equation) fully describes a single projection ray, since it exactly describes the list of pixels that are hit by that ray. Each ray will only meet a very small percentage of the pixels, since it crosses the reconstruction area in a straight line. Hence, the matrix \({\bf A}\) is very sparse (a large percentage of its elements is zero). If the reconstruction grid is square (say, 512×512 pixels, where I’m assuming 2D reconstruction), and you take 360 projections with a detector of 725 (about \(\sqrt{2}\times 512\), why?) pixels long, then \({\bf A}\) will be 261000×262144 (more than 68 billion elements!). It helps that it is sparse, but, in practice, the matrix is almost never computed explicitly. The operation of the matrix is simply recomputed each time it is needed. Writing a computer program that creates a matrix like this in practice is not really difficult, but will take you at least a few hours, and make you dream of sines and cosines the following night.

Algebraic Reconstruction

For the algebraic reconstruction itself, there are lots of possibilities. For starters, you can treat it as a general linear system of equations, and solve it using one of the many available techniques. In fact, some of the techniques that were developed specifically for tomography were later discovered to be reinventions of existing methods.

Another way of solving \({\bf Ax}={\bf b}\) is the Moore-Penrose pseudoinverse \({\bf A}^+\). The solution is then simply \({\bf x}^+={\bf A}^+{\bf b}\). But there are several problems with this. First, the matrix is almost always simply too large to compute the pseudoinverse explicitly. Second, even if that is be possible, \({\bf A}\) is often almost singular, which makes it numerically difficult to compute \({\bf A}^+\).

In practice, the algebraic reconstruction problem is typically attacked using iterative techniques. These start from a certain initial guess of the solution, and then repeatedly apply an update step, till the solution is “good enough” according to some criterion. Most, but not all, of these iterative methods (mathematically) converge to a solution. Some methods are based on heuristic (experience-based) arguments that have no mathematically guaranteed convergence, but are still very useful in practice. Examples of iterative algorithms are SIRT and PDART.

This is the last article on tomography, for the time being. These four articles have, of course, only scratched the surface of this broad and fascinating subject. But, for now, I leave it at this…