Axial Tomography by filtered least squares

Kenneth C. Holmes

Max Planck Institute for Medical Research, 69120 Heidelberg, Germany

e-mail – holmes@mpimf-heidelberg.mpg.de

and

Rasmus R. Schroeder

CellNetworks, BioQuant, University of Heidelberg

Im Neuenheimer Feld 267, 69120 Heidelberg, Germany

e-mail  – rasmus.schroeder@bioquant.uni-heidelberg.de

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 A₀ has components

h1=1, h2=0, h3=0, h4=0 etc.

and for B₀

h₁ = 0 h₅ = 1 h₉ = 0

h₂ = 1 h₆ = 1 h₁₀ = 0

h₃ = 1 h₇ = 1 h₁₁ = 0

h₄ = 1 h₈ = 1 h₁₂ = 0 etc.

Similar results are obtained by inspection for each of the other points C,D,E,F etc.

Fig. 1

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

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 x = b (1)

where b is the list of all observables i.e. all the projected densities A0,B0, ... A1,B1... A2,B2....A12, B12 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.

The Solution

The classical solution to the equations

H x= b+ ε (1a)

which minimizes the sum of squares of errors e^Te (least squares solution) is

H^TH x = H^Tb (2)

H^TH is the normal or Hessian matrix, which from its definition is square symmetric. If H^TH is non-singular it can be inverted to give the solution

x = (H^TH)^-1 H^T b (3)

However, in our case b is shorter than x (we do not have enough observations to determine the density x). Because of the form of H, the vector H^Tb is longer than b. Normally, in least squares there is a reduction of dimensionality on pre-multiplying by H. However, in the present case since the number of unknown exceeds the number of observations the dimensionality of the normal matrix is actually increased. Equation 2 is ill determined and cannot be solved by normal methods.

We proceed by expressing H^TH in diagonal form. Since H^TH is square and positive semi-definite one can find a transformation

H^TH = V^TΛV (4)

where Λ is a diagonal matrix of the eigenvalues λi. (the eigenvalues of H^TH). V is a matrix of which the columns are the orthonormal eigenvectors of H^TH. Premutiplication by V rotates into a representation where the eigenvectors become the principle axes of the problem. Substituting (4) in (3) we have

V^TΛVx = H^Tb (5)

or ΛVx= VH^Tb (6)

We are now in a representation where the principle axes of the Hessian or normal matrix are the base and all equations have been separated. Now we may solve (6) element by element by dividing by λi

(Vx)_i = 1/λ_i (V H^Tb)_i (7)

We can simply stop solving for elements of Vx when λ gets small compared with the errors in the observations. Hence we determine some (but not all) values of Vx. The rest we set to zero. We call this truncated vector (Vx)’

Rotating (transforming) back gives

x = V^T (Vx)’ (8)

that is the required solution.

We actually start out with the back projection

H^Tb, the right hand side of equation 2, is a vector derived from the observations that is central feature of the least squares method. Noteworthy, in the present context, H^Tb assumes a special significance; it is in fact the back projection! This comes about because H is a projector matrix.

The Eigenfunctions

Each of the eigenvectors of H^TH (one of the the columns of V) may be mapped onto the lattice since each component of the vector is associated with a lattice point. When we do this we generate two-dimensional functions in the circular space we have chosen for the solution. These functions are (by definition) orthonormal. Because we have chosen a circular boundary, they look rather like cylinder functions. In fact they have a symmetry characteristic of the geometry of the problem. Some examples from a case considered below are shown in Fig.3

Fig 3

A transformation and back transformation yields the answer

The significance of equations 6, 7 and is as follows: equation 6 is a transformation of the back projection using a set of orthonormal functions that are the eigenfunctions of the outer product of the projector matrix with itself. The result of this tranformation is a list of numbers (with associated errors). To proceed we divide each of these numbers by the appropriate eigenfunction arranged from large to small in order to determine how much of each eigenfunction contribute to the final answer (7). At some point the eigenfunction becomes as small as the errors and we have to stop, otherwise we amplify up noise – there is no point in going on – a truncation. This yields a list of numbers telling us how much of each eigenfunction to put into the answer. The result is calculated by a back transformation (8). The effect of truncation is to limit the resolution in the answer. It transpires that the eigenfunctions with the longest spatial frequencies have the largest eigenvalues so that the resolution goes up with more eigenvalues. Since the base we have used is orthonormal and final we would not get another answer if we considered a larger subspace – just more noise.

A Donkey

Fig 4

To test the theory, we took a donkey as a trial object. This object is demanding in a number of respects. It has no symmetry, It has sharp graduations of density. It has long legs! For first trial the projected density was calculated every 3º between 0 and 180º (the symmetry of the projection precludes the addition of any new information in the second half circle). Projections at 0º and 30º are shown in Fig 4. The 61 projections were used to reconstruct the donkey using the theory given above. The diameter of the donkey was about 150 pixels, so we set up a circle of radius 75 pixels. This produced 17645 lattice points (π.75² = 17671). Thus the order of the normal matrix is 17645. Calculation shows that the effective rank is 5000-6000 (the cut off is not sharp). Using a Fortran program running on a Silicon Graphics Origin 3000 and four processors it took about 1/2h. to set up the normal matrix and another 3 h. to extract the 4000 most significant eigenvectors using the Lapack routine SSYEVX. It is noteworthy that this considerable computing load has only to be undertaken for a new geometry. If the geometry remains unaltered (same number of projections in the same angular range) then a reconstruction can be calculated from stored matrices in a fraction of a second. In Fig 5 we show the results of this reconstruction. Shown are; the original object; the raw back projection, the r-weighted back-projection, the reconstructed image using 2000 largest eigenfunctions; the reconstructed image using 4000 eigenfunctions.

Fig 5

The use of a mask

We are not obliged to use a round area to define the boundaries of the lattice. If we have some knowledge of the shape of the object to be reconstructed this can be used as a mask to define the area within which we seek the density of the object. The use of a mask will reduce the size of the normal matrix, which has two desirable effects. The calculation of the diagonalized matrix is much faster; and the resolution obtainable from the given data is improved since we are seeking to solve for a smaller number of variables. The use of a mask, which fits the object, is rather like the use of solvent flattening in protein crystallography to improve the resolution of electron density maps. It the following example the projections have been made every 6º, giving 30 independent projections. Fig.6 shows the back projection and r-weighted back projection. Fig 7, shows the effects of the mask.

This method has been successfully used to reconstruct actin decorated with myosin cross-bridges at 13Å resolution [5].

Fig 6 Back projection and r-weighted back projection for 30 projections at 6º intervals

Fig 7. (see next page) Projections at 6º intervals. The effects on resolution of a mask compared with a circle for a given number of eigenfunction. The eigenvalues fall off more quickly for the mask than for the circle. Thus the circle eigenvalues are down to 1% of the origin value at 1300 eigenvalues and 0.1% at 2300 eigenvalues. The corresponding numbers for the mask are 1000 and 1800.

Fig. 7

References

1. Klug, A. and R.A. Crowther, Three-dimensional image reconstruction from the viewpoint of information theory. Nature, 1972. 238: p. 435-440.

2. Gilbert, P. Ph.D. thesis 1972 University of Cambridge

3. Crowther, R.A., D.J. DeRosier, and A. Klug, The reconstruction of a three-dimensional structure from projections and its application to electron microscopy. Proceedings of the Royal Society of London Series A Mathematical & Physical Sciences, 1972. 317: p. 319-40.

4. Diamond, R., Filtering in the method of least squares. Acta Cryst., 1969. B25: p. 805-806.

5. Holmes, K. C., Angert, I., Kull, F. J., Jahn, W. and Schroder, R. R. 2003. Electron cryo-microscopy shows how strong binding of myosin to actin releases nucleotide. Nature 425: p. 423-7.