International
Tables for
Crystallography
Volume B
Reciprocal space
Edited by U. Shmueli

International Tables for Crystallography (2010). Vol. B, ch. 2.5, pp. 369-371   | 1 | 2 |

Section 2.5.6.4. The algebraic and iterative methods

B. K. Vainshteinc and P. A. Penczekg

2.5.6.4. The algebraic and iterative methods

| top | pdf |

The algebraic methods have been derived based on the observation that when the projection equation (2.5.6.5)[link] is discretized, it forms a set of linear equations. Thus, pixels from all available N projections are placed (in an arbitrary order) in a vector [\varphi _{2jk}^{\,n} , n = 1, \ldots ,N \to {\bf{f}}] and the voxels of the 3D object in a vector [\varphi _{3jkl} \to {\bf{g}}] (in an order derived from the order of f by algebraic relations). Note we left the exact sizes of f and g undetermined, as the major advantage of algebraic methods is that we can include in the reconstruction only pixels located within an arbitrary support in two dimensions and this support can be different for each projection; similarly, the support of the object in three dimensions can be arbitrary. Thus, the number of elements in f and g are at most K2N and K3, respectively, with K2 being the number of pixels within chosen support. In the algebraic formulation the operation of projection is defined by the projection matrix P whose elements [p_{\eta \xi } ] are the interpolation weights. Their values are determined by the interpolation scheme used, but for the bi- and trilinear interpolations [0 \le p_{\eta \xi } \le 1]. The algebraic version of (2.5.6.5)[link] is[{\bf{f}} = {\bf{Pg}}.\eqno(2.5.6.19)]Matrix P is rectangular and since in single-particle reconstruction the number of projections exceeds the linear size of the object ([N \gg K]) the system of equations is overdetermined. It can be solved in a least-squares sense:[{\tilde{\bf g}} = ( {{\bf{P}}^T {\bf{P}}} )^{ - 1} {\bf{P}}^T {\bf{f}},\eqno(2.5.6.20)]which yields a unique structure [{\tilde{\bf g}}] that corresponds to the minimum of [| {{\bf{Pg}} - {\bf{f}}} |^2 ]. As in the case of direct inversion in Fourier space (2.5.6.17)[link], the approach (2.5.6.20)[link] is impractical because of the very large size of the projection matrix. Indeed, the size of P is K2N × K3, which for a modest number of projections N = 10 000 and image size K2 = 642 = 4096 yields a matrix of size ~(4 × 107) × (2 × 106)! Nevertheless, in some cases the projections of the structure are orthoaxial, the full 3D reconstruction reduces to a series of independent 2D reconstructions, and it becomes possible to solve (2.5.6.19)[link] by using the singular value decomposition (SVD) of the matrix [{\bf{P}}]:[{\bf{P}} = {\bf U}\boldSigma {\bf V}^T,\eqno(2.5.6.21)]where [{\bf{U}}] and [{\bf{V}}] contain eigenvectors of matrices [{\bf{PP}}^T ] and [{\bf{P}}^T {\bf{P}}], respectively, and [\boldSigma] is a diagonal matrix containing the first r nonnegative square roots of the eigenvalues of [{\bf{PP}}^T ] (Golub & Van Loan, 1996[link]). The solution is given by[{\bf{g}} = {\bf{V}}\boldSigma _l^{ - 1} {\bf{U}}^T {\bf{f}},\eqno(2.5.6.22)]where [\boldSigma_l^{ - 1} ] is the inverse of the matrix containing the [l \le r] largest elements of [\boldSigma] with the remaining values set to zero. By selecting l appropriately we achieve a measure of regularization. The advantage of the approach is that for a given geometry the decomposition has to be calculated once; thus, the method becomes very efficient if the reconstruction has to be performed repeatedly for the same distribution of projections or if additional symmetries, such as helical, are taken into account (Holmes et al., 2003[link]).

In the general 3D case, a least-squares solution can be found using one of the iterative approaches that take advantage of the fact that the projection matrix is sparse. Indeed, if bilinear interpolation is used, a row of [{\bf{P}}] will contain only about 4N nonzero elements (out of the total ~K2N). The main idea is that the matrix [{\bf{P}}] is not explicitly calculated or stored; instead, its elements are calculated repeatedly during iterations as needed. In the simultaneous iterative reconstruction technique (SIRT) we find the minimum of[L( {\bf{g}} ) = | {{\bf{Pg}} - {\bf{f}}} |^2\eqno(2.5.6.23)]by selecting the initial 3D structure [{\bf{g}}^0 ] (usually set to zero) and by iteratively updating its current approximation [{\bf{g}}^{i + 1} ] using the gradient of the objective function (2.5.6.23)[link] [\nabla L( {\bf{g}} )],[{\bf{g}}^{i + 1} = {\bf{g}}^i - \lambda ^i {\bf{P}}^T ( {{\bf{Pg}}^i - {\bf{f}}} ) = {\bf{g}}^i - \lambda ^i ( {{\bf{P}}^T {\bf{Pg}}^i - {\bf{P}}^T {\bf{f}}} ).\eqno(2.5.6.24)]Setting the relaxation parameter [\lambda ^i = \lambda = {\rm constant}] yields Richardson's algorithm (Gilbert, 1972[link]). In single-particle reconstruction applications convergence (i.e., a solution that is not dominated by the noise) is usually reached in ~100 iterations.

SIRT is extensively used in single-particle reconstruction (Frank et al., 1996[link]) because it yields superior results over a wide range of experimental conditions (Penczek et al., 1992[link]) and, when presented with angular gaps in the distribution of projections, produces the least disturbing artifacts. Note that strictly speaking, filtered backprojection and direct Fourier inversion algorithms are not applicable to data that do not cover the Fourier space fully. In addition, SIRT offers considerable flexibility in EM applications. First, it is possible to accelerate the conver­gence by adjusting relaxation parameters: setting [\lambda ^i =] [ \arg \min \nolimits_{\lambda \ge 0} L( {g^i - \lambda \nabla L( {g^i } )} )] results in a steepest decent algorithm. Secondly, even faster convergence (in ~10 iterations) is achieved by solving (2.5.6.23)[link] using the conjugate-gradient method, but this requires addition of a regularizing term to (2.5.6.23)[link] in order to prevent excessive enhancement of noise. Such a term has the form [| {{\bf{Bg}}} |^2 ], where matrix B is a discrete approximation of a Laplacian or higher-order derivatives. Thirdly, it is possible to take into account the underfocus settings of the microscope by including the contrast transfer function (CTF) of the microscope in (2.5.6.23)[link] and by solving the problem for the structure g that in effect will be corrected for the CTF,[L( {\bf{g}} ) = ( {1 - \eta } )\textstyle\sum\limits_\kappa {| {{\bf{S}}_\kappa {\bf{P}}_\kappa {\bf{g}} - {\bf{f}}_\kappa } |^2 } + \eta | {{\bf{Bg}}} |^2 .\eqno(2.5.6.25)]Here we assumed that the projection data were grouped according to the defocus settings indexed by [\kappa ] and we introduced a Lagrange multiplier [\eta ] whose value determines the smooth­ness of the solution g (Zhu et al., 1997[link]). [{\bf{S}}_\kappa ] represents the space-invariant point-spread function of the microscope for the κth defocus group; thus [{\bf{S}}_\kappa ] has a block-Toeplitz structure (Biemond et al., 1990[link]) and [{\bf{S}}_\kappa = {\bf{S}}_\kappa ^T ]. If we assume that there is no astigmatism (which means that the point-spread function is rotationally invariant), [{\bf{S}}_\kappa ] can be either applied to the projection of the structure or to the 3D structure itself. With these assumptions, (2.5.6.25)[link] can be solved using the iterative scheme (Penczek et al., 1997[link])[\eqalignno{{\bf g}^{i + 1} &= {\bf g}^i - \lambda\big\{ ( 1 - \eta )\big[ \textstyle\sum\limits_\kappa ( {\bf S}_\kappa ^T {\bf P}_\kappa ^T {\bf P}_\kappa {\bf S}_\kappa {\bf g}^i ) - \textstyle\sum\limits_\kappa ( {\bf S}_\kappa ^T {\bf P}_\kappa ^T {\bf f}_\kappa ) \big] &\cr &\quad+ \eta {\bf B}^T {\bf Bg}^i \big\}.&(2.5.6.26)}]The second sum can be precalculated and stored as a 3D volume in the computer memory. Thus, the input projections have to be read only once and are never accessed again during the course of the iterations, which eliminates the need to store them in the memory [as is also the case in (2.5.6.24)[link]]. In addition, the product BTB is the 3D Laplacian, which can be applied efficiently without actually creating the matrix BTB. Finally, because of the large size of the matrix [{\bf{S}}_\kappa ] it is more convenient to apply it in Fourier space and to modify the Fourier transform of the volume instead of using matrix multiplication or real-space convolution. Algorithms (2.5.6.24)[link]–(2.5.6.26)[link][link] are implemented in the SPIDER package (Frank et al., 1996[link]) and (2.5.6.24)[link] in the SPARX package (Hohn et al., 2007[link]).

The algebraic reconstruction technique (ART) predates SIRT; in the context of tomographic reconstructions it was proposed by Gordon and co-workers (Gordon et al., 1970[link]) and later it was recognized as a version of Kaczmarz's method for iteratively solving (2.5.6.23)[link] (Kaczmarz, 1993[link]). We write (2.5.6.19)[link] as a set of systems of equations, each relating single pixels fn, [n = 1, \ldots ,NK^2 ] in projections with the 3D structure,[f_n = {\bf{p}}_n^T {\bf{g}},\quad n = 1, \ldots ,NK^2.\eqno(2.5.6.27)]Note (2.5.6.27)[link] and (2.5.6.19)[link] are equivalent, as [{\bf{f}} = [ {f_1 f_2 \ldots f_N } ]^T ] and [{\bf{P}} = [ {{\bf{p}}_1 {\bf{p}}_2 \ldots {\bf{p}}_N } ]^T ]. With this notation and a relaxation parameter [0 \,\lt\, \mu \,\lt\, 2], ART comprises the following steps:

  • (1) set i = 0 and the initial structure [{\bf{g}}^0 ];

  • (2) for [n = 1, \ldots ,NK^2 ], set[{\bf{g}}^{iNK^2 + n} = {\bf{g}}^{iNK^2 + n - 1} - \mu ( {{\bf{p}}_{{n}}^{{T}} {\bf{g}}^{iNK^2 + n - 1} - f_n } )({{{\bf{p}}_n }/{| {{\bf{p}}_n } |}})\semi\eqno(2.5.6.28)]

  • (3) set i = i + 1; go to step (2)[link].

Although the mathematical forms of update equations in SIRT (2.5.6.24)[link] and in ART (2.5.6.28)[link] are very similar, there are profound differences between them. In SIRT, all voxels in the structure are corrected simultaneously after projections and backprojections of the current update of the structure are calculated. In ART, the projection/backprojection step (2.5.6.28)[link] involves only correction with respect to an individual pixel in a single projection immediately followed by the update of the structure. This results in a much faster convergence of ART as compared to SIRT. Further acceleration can be achieved by selecting the order in which pixels enter the correction step (2)[link] in (2.5.6.28)[link]. It was observed that if a pixel is selected such that its projection direction is perpendicular to the projection direction of the previous pixel, the convergence is achieved faster (Hamaker & Solmon, 1978[link]; Herman & Meyer, 1993[link]). Interestingly, a random order works almost equally well (Natterer & Wübbeling, 2001[link]).

In single-particle reconstruction, ART has been introduced in the form of `ART with blobs' (Marabini et al., 1998[link]) and is available in the Xmipp package (Sorzano et al., 2004[link]). In this implementation, the reconstruction structure is represented by a linear combination of spherically symmetric, smooth, spatially limited basis functions, such as Kaiser–Bessel window functions (Lewitt, 1990[link], 1992[link]; Matej & Lewitt, 1996[link]). Introduction of blobs significantly reduces the number of iterations necessary to reach an acceptable solution (Marabini et al., 1998[link]).

The major advantage of iterative reconstruction methods is the ability to take advantage of a priori knowledge, i.e., any information about the protein structure that was not initially included in the data processing, and introduce it into the reconstruction process in the form of constraints. Examples of such constraints include similarity to the experimental (measured) data, positivity of the protein mass density (only valid in conjunction with the CTF correction), bounded spatial support etc. Formally, the process of enforcing selected constraints is best described in the framework of the projections onto convex sets (POCS) theory (Sezan & Stark, 1982[link]; Youla & Webb, 1982[link]; Sezan, 1992[link]) introduced into EM by Carazo and co-workers (Carazo & Carrascosa, 1986[link], 1987[link]; Carazo, 1992[link]).

References

Biemond, J., Lagendijk, R. L. & Mersereau, R. M. (1990). Iterative methods for image deblurring. Proc. IEEE, 78, 856–883.
Carazo, J. M. (1992). The fidelity of 3D reconstruction from incomplete data and the use of restoration methods. In Electron Tomography, edited by J. Frank, pp. 117–166. New York: Plenum.
Carazo, J. M. & Carrascosa, J. L. (1986). Information recovery in missing angular data cases: an approach by the convex projections method in three dimensions. J. Microsc. 145, 23–43.
Carazo, J. M. & Carrascosa, J. L. (1987). Restoration of direct Fourier three-dimensional reconstructions of crystalline specimens by the method of convex projections. J. Microsc. 145, 159–177.
Frank, J., Radermacher, M., Penczek, P., Zhu, J., Li, Y., Ladjadj, M. & Leith, A. (1996). SPIDER and WEB: processing and visualization of images in 3D electron microscopy and related fields. J. Struct. Biol. 116, 190–199.
Gilbert, P. F. C. (1972). Iterative methods for the three-dimensional reconstruction of an object from projections. J. Theor. Biol. 36, 105–117.
Golub, G. H. & Van Loan, C. F. (1996). Matrix Computations. Baltimore: Johns Hopkins University Press.
Gordon, R., Bender, R. & Herman, G. T. (1970). Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography. J. Theor. Biol. 29, 471–481.
Hamaker, C. & Solmon, D. C. (1978). Angles between null spaces of X-rays. J. Math. Anal. Appl. 62, 1–23.
Herman, G. T. & Meyer, L. B. (1993). Algebraic reconstruction techniques can be made computationally efficient. IEEE Trans. Med. Imaging, 12, 600–609.
Hohn, M., Tang, G., Goodyear, G., Baldwin, P. R., Huang, Z., Penczek, P. A., Yang, C., Glaeser, R. M., Adams, P. D. & Ludtke, S. J. (2007). SPARX, a new environment for cryo-EM image processing. J. Struct. Biol. 157, 47–55.
Holmes, K. C., Angert, I., Kull, F. J., Jahn, W. & Schroder, R. R. (2003). Electron cryo-microscopy shows how strong binding of myosin to actin releases nucleotide. Nature (London), 425, 423–427.
Kaczmarz, S. (1993). Approximate solutions of systems of linear equations (Reprint of Kaczmarz, S., Angenäherte Auflösung von Systemen linearer Gleichunger, Bulletin International de l'Academie Polonaise des Sciences. Lett A, 355–357, 1937). Int. J. Control, 57, 1269–1271.
Lewitt, R. M. (1990). Multidimensional digital image representations using generalized Kaiser–Bessel window functions. J. Opt. Soc. Am. A, 7, 1834–1846.
Lewitt, R. M. (1992). Alternatives to voxels for image representation in iterative reconstruction algorithms. Phys. Med. Biol. 37, 705–716.
Marabini, R., Herman, G. T. & Carazo, J. M. (1998). 3D reconstruction in electron microscopy using ART with smooth spherically symmetric volume elements (blobs). Ultramicroscopy, 72, 53–65.
Matej, S. & Lewitt, R. M. (1996). Practical considerations for 3-D image reconstruction using spherically symmetric volume elements. IEEE Trans. Med. Imaging, 15, 68–78.
Natterer, F. & Wübbeling, F. (2001). Mathematical Methods in Image Reconstruction. Philadelphia: SIAM.
Penczek, P., Radermacher, M. & Frank, J. (1992). Three-dimensional reconstruction of single particles embedded in ice. Ultramicroscopy, 40, 33–53.
Penczek, P. A., Zhu, J., Schröder, R. & Frank, J. (1997). Three-dimensional reconstruction with contrast transfer function compensation from defocus series. Scan. Microsc. Suppl. 11, 1–10.
Sezan, M. I. (1992). An overview of convex projections theory and its application to image recovery problems. Ultramicroscopy, 40, 55–67.
Sezan, M. I. & Stark, H. (1982). Image restoration by the method of convex projections. II. Applications and numerical results. IEEE Trans. Med. Imaging, 1, 95–101.
Sorzano, C. O. S., Marabini, R., Velazquez-Muriel, J., Bilbao-Castro, J. R., Scheres, S. H. W., Carazo, J. M. & Pascual-Montano, A. (2004). XMIPP: a new generation of an open-source image processing package for electron microscopy. J. Struct. Biol. 148, 194–204.
Youla, D. C. & Webb, H. (1982). Image restoration by the method of convex projections. 1. Theory. IEEE Trans. Med. Imaging, 1, 81–94.
Zhu, J., Penczek, P. A., Schröder, R. & Frank, J. (1997). Three-dimensional reconstruction with contrast transfer function correction from energy-filtered cryoelectron micrographs: procedure and application to the 70S Escherichia coli ribosome. J. Struct. Biol. 118, 197–219.








































to end of page
to top of page