# Axial Tomography by filtered least squares

#### Introduction

We address the problem of reconstructing the 2-dimensional density of a plane object from a number of projections of the density in the plane of the object. The one-dimensional projections are generated by rotating the object around an axis at right angles to the plane and projecting down a constant direction in the plane of the object. Thus the direction of projection is defined to be at right angles to the rotation axis. i.e. the geometry of the single axis projection. Klug and DeRosier [1] first solved this problem by showing how the two-dimensional Fourier transform of the required density could be built up from transforms of the projected density. This is the basis of the method of reciprocal-space analysis and Fourier synthesis of the density. Another frequently used method makes use of back projection, where each projection is used to generate a set of lines along the direction of projection, each line weighted by the projected density. The sum over all such lines for each point in the object gives an estimate of the reconstructed density. Back projection yields an image in which the high frequency components are attenuated. A way to compensate for this is to compute the Fourier transform of the back projection (with its origin on the axis of rotation) and to sharpen the Fourier transform of the density computed by back projection by a weighting function that is proportional to the distance from the origin in Fourier space (Gilbert [2]). Thus the center of the transformed back projection receives zero weight, and at some resolution there is a sharp cut off.

Crowther DeRosier and Klug [3] examined the possibility of using a real space method of reconstructing the required density based on filtered least squares. One defines the sought-after density on a grid and various projections of this grid are available as observations. This problem can then be solved by least squares. Moreover, an algebraic formulation allows an analysis of the information content of the available projections. Typically there will be too few projection equations to determine the problem at the resolution of the grid. The rank of the normal matrix will be less than the order and therefore the normal matrix has no unique inverse.

The problem actually goes deeper: the rotation axis introduces a singularity into the matrix that cannot be removed by adding more views. This singularity has bedeviled all algebraic formulations of the axial tomography problem. However, as was suggested by Crowther DeRosier and Klug (loc. cit.) one can get round these problems by diagonalizing the normal matrix and by inverting the normal matrix in the subspace of significant eigenvalues. This method of inverting ill-defined matrices was described by Diamond [4]. The truncation of the inversion limits the resolution in the final density. However, this filtered least squares was never used as a method of reconstruction because at the time the computing load would have been prohibitive. However, computer power in no longer limiting and the method is now applicable. In the following we show how the density of the object may be recovered from the projected density using the method of least squares. The method extracts the maximum possible amount of information about the reconstructed density in an unbiased way. Moreover, the crude shape of the object can be used to define the shape of the domain in which the solution is sought. We refer to this as a mask. The use of a mask (which is rather like solvent flattening on crystallography) yields more accurate reconstructions. It is particularly valuable for asymmetric objects. Furthermore, it can significantly reduce the order and rank of the normal matrix, which is computationally very advantageous.

Finally, it turns out that the right hand side of the normal equations is the back projection. Thus the theory demonstrates that an analytical procedure exists for extracting the three dimensional density from the back projection. Indeed, the correct filter for turning the back projection into the reconstructed object is in fact the inverted normal matrix!

The method we have developed has close parallels with the weighted back projection. However, in place of the sine-cosine transform we operate on the back projection with an orthogonal matrix of eigenvectors derived from the normal matrix. The nature of the eigenvectors (and the derived eigenfunctions) depends on the shape of the mask and the number and spacing of the projections. For a round mask they look like cylinder functions. This operation is analogous to a Fourier transformation. In the transformed space the coefficients are divided by the appropriate calculated eigenvalue: the eigenvalues are large for the low orders and fall off monotonically. Thus the transformed density is subjected to a sharpening analogous to r weighting (Gilbert [2]). Somewhere, because the eigenvalues get very small, we have to truncate. The point depends on the noise in the data. Thus we see that noise limits resolution. Finally, we operate on the sharpened transform with the inverse orthogonal matrix of eigenvectors, which is analogous to a back Fourier transform. Compared with the weighted back projection method the least-squares method uses custom-made orthogonal functions for computing the transform of the back projection. These functions are determined by and are characteristic of the geometry of the problem. Moreover, the weighting function does not go to zero at the origin of the transform and is also appropriate for the geometry. Thus the method deals quantitatively with the relationship between low and high order terms. The method yields the best possible solution for the given geometry without any ad hoc assumptions and moreover, provides a theoretical basis for analyzing the weighted back projection method.

#### Setting up the problem

We represent the unknown density in the selected plane as values on a lattice within some chosen radius (Fig 1a). The use of a circular boundary is arbitrary. However, round is the best choice in the absence of prior information about the object. The lattice points are tagged with serial numbers so that they can be written down as a long list of numbers (a vector x). For example, in the first observation (rotation zero) the values at the points 2-8 add up to the projected density B0. The points 9-19 sum to the projected density C0 etc. We can write the projector for each point A0, B0 etc. as a vector h, which spans all points. In the present case by inspection the projector for A0 has components

h_{1}=1, h_{2}=0, h_{3}=0, h_{4}=0 etc.

and for B_{0}

h_{1}=0 h_{5}=1 h_{9}=0

h_{2}=1 h_{6}=1 h_{10}=0

h_{3}=1 h_{7}=1 h_{11}=0

h_{4}=1 h_{8}=1 h_{12}=0 etc.

**Fig 1**: In the second measurement (plane rotated for example by 360/13º) the situation is not quite so straightforward

**Fig 1**: In the second measurement (plane rotated for example by 360/13º) the situation is not quite so straightforward

© Kenneth C. Holmes

Similar results are obtained by inspection for each of the other points C,D,E,F etc. In the second measurement (plane rotated for example by 360/13º) the situation is not quite so straightforward (fig 1b). Again the integration is carried out by sampling the density at equal intervals along the projecting lines but the integration points (shown as small circles) no longer coincide with the lattice points. We can no longer write down the projector by inspection. An interpolation scheme is needed. We have used Gaussian interpolation. Each lattice point in the plane is convoluted with a gauss function of suitable width (Fig 2).

**Fig 2**: An interpolation scheme is needed. We have used Gaussian interpolation. Each lattice point in the plane is convoluted with a gauss function of suitable width.

© Kenneth C. Holmes

As before, integration is carried out at equal distances along projector lines but at each integration point the contributions of all nearby lattice points are taken into account (weighted according to the value of the appropriate Gaussian at the integration point). Thus a number of lattice points will contribute to each integration point. By extending this method to all projected points for each rotational position we finally obtain observational equations relating all the unknown density points with the known projections, which we can write as

H_{4}=**b** (1)

where **b** is the list of all observables i.e. all the projected densities A_{0}, B_{0}, ... A_{1}, B_{1}, ... A_{2}, B_{2} .... A_{12}, B_{1} in order and **x** is a vector of the unknown weights of the Gaussians sitting on each of the defined lattice points which we wish to determine. H is a projector matrix derived by assembling all the available vectors **h** arranged in rows. Each row is as long as the list of unknown density points and consists mainly of zero's. The number of rows is equal to the number of observations, so that H is very wide and not very high! We wish to determine **x**.