International
Tables for Crystallography Volume B Reciprocal space Edited by U. Shmueli © International Union of Crystallography 2010 
International Tables for Crystallography (2010). Vol. B, ch. 2.5, pp. 366375

2.5.6. Threedimensional reconstruction^{5}
In electron microscopy (EM) and singleparticle reconstruction, threedimensional (3D) reconstruction methods are used for studying biological structures; that is, symmetric or asymmetric associations of biomacromolecules (muscles, spherical and rodlike viruses, bacteriophages, individual proteins and ribosomes) (Frank, 2006). The electron microscope is used to obtain parallelbeam twodimensional (2D) projections (ray transform) of frozen hydrated 3D macromolecules suspended in random orientations (Fig. 2.5.6.1). The function is the 2D projection of the 3D molecular electron distribution . One can also consider onedimensional (1D) projections of multidimensional distributions; the set of these projections is called a Radon transform. For 2D distributions, ray and Radon transforms differ only in the notation. For (d > 2)dimensional distributions the two are different: in a ray transform the integrals are calculated over straight lines and yield (d − 1)dimensional projections, while in Radon transforms the integrals are calculated over (d − 1)dimensional hyperplanes and yield 1D projections. In electron microscopy, Radon transforms are not directly measurable, but can be formed computationally and used in intermediate steps of the 3D reconstruction or in alignment procedures (Radermacher, 1994).

A threedimensional object and its twodimensional projection . The electron beam penetrates the specimen in the direction of the z axis. 
Within the linear weakphaseobject approximation of the imageformation process in the microscope [see equation (2.5.2.42) in Section 2.5.2 of this chapter], 2D projections represent line integrals of the potential of the particle under examination convoluted with the pointspread function of the microscope, s, so, using (2.5.2.43),Sincewe haveIf we omit constant terms, we obtain
In this section, we will assume that all images were collected using the same defocus setting, so the pointspread function s is constant and does not depend on the projection direction . Thus, we will concern ourselves with the inversion of the projection problemThe projection direction is defined by a unit vector and it is formed on the plane x perpendicular to . The set of various projections may be assigned by a discrete or continuous set of points on a unit sphere (Fig. 2.5.6.2).

The projection sphere and projection of along τ onto the plane . The case represents orthoaxial projection. Points indicate an arbitrary distribution of projection directions τ. 
In Fourier space, the relation between an object and its projection is referred to as the central section theorem: the Fourier transformation of projection of a 3D object is the central (i.e., passing though the origin of reciprocal space) 2D plane cross section of a 3D transform perpendicular to the projection vector (Bracewell, 1956; DeRosier & Klug, 1968; Crowther, Amos et al., 1970). In Cartesian coordinates, a 3D Fourier transform isThe transform of projection along z isIn the general case of projecting along the vector , the central section theorem isFrom this theorem it follows that the inversion of the 3D ray transform is possible if there is a continuous set of projections corresponding to the motion of the vector over any continuous line connecting the opposite points on the unit sphere (Fig. 2.5.6.2) (Orlov's condition: Orlov, 1976). This result is evidenced by the fact that in this case the cross sections that are perpendicular to in Fourier space continuously cover the whole Fourier space, i.e., they yield and thereby determine .
In singleparticle reconstruction, imaged objects are randomly and nonuniformly oriented on the substrate at different angles (Frank, 2006) and the distribution of their orientations is beyond our control; therefore, the practical impact of Orlov's condition is limited. In fact, it is more useful to determine a posteriori, i.e., after the 3D reconstruction of the macromolecule is computed, how well the Fourier space was covered. This can be done by calculating the distribution of the 3D spectral signaltonoise ratio (Penczek, 2002).
In the general case of the 3D reconstruction of from projections , the projection vector occupies arbitrary positions on the projection sphere (Fig. 2.5.6.2). First, let us consider the case of a 2D function and its ray transform . We introduce an operation of backprojection b, which is stretching along each 1D projection (Fig. 2.5.6.3). When the result is integrated over the full angular range of projections , we obtain the projection synthesis defined asHowever, the backprojection operator is not the inverse of a 2D ray transform, as the resulting image b is blurred by the pointspread function (Vainshtein, 1971):By noting that the Fourier transform of is and by using the convolution theorem , we obtain the `backprojectionfiltering' inversion formula:The more commonly used `filteredbackprojection' inversion is based on the 2D version of the central section theorem (2.5.6.8):where . With this in mind, can be related to its ray transform by evaluating the Fourier transform of in polar coordinates:In three dimensions, the backprojection stretches each 2D projection along the projection direction . A 3D synthesis is the integral over the hemisphere (Fig. 2.5.6.2)Thus, in three dimensions the image b obtained using the backprojection operator is blurred by the pointspread function (Vainshtein, 1971). It is possible to derive inversion formulae analogous to (2.5.6.11) and (2.5.6.13).

(a) Formation of a backprojection function; (b) projection synthesis (2.5.6.9) is a superposition of these functions. 
The inversion formulae demonstrate that it is possible to invert the ray transform for continuous functions and for a uniform distribution of projections. In electron microscopy, the projections are never distributed uniformly in three dimensions. Indeed, a uniform distribution is not even desirable, as only certain subsets of projection directions are required for the successful inversion of a 3D ray transform, as follows from the central section theorem (2.5.6.8). In practice, we always deal with sampled data and with discrete, random and nonuniform distributions of projection directions. Therefore, the inversion formulae can be considered only as a starting point for the development of the numerical (and practical!) reconstruction algorithms. According to (2.5.6.10) and (2.5.6.14), a simple backprojection step results in reconstruction that corresponds to a convolution of the original function with a pointspread function that depends only on the distribution of projections, but not on the structure itself. Taking into account the linearity of the backprojection operation, one has to conclude that for any practically encountered distribution of projections it should be possible to derive the corresponding pointspread function and then, using either deconvolution or Fourier filtration (with a `weighting function'), arrive at a good approximation of the structure. This observation forms the basis of the weighted backprojection algorithm (Section 2.5.6.5). Similarly, the central section theorem gives rise to direct Fourier inversion algorithms (Section 2.5.6.6). Nevertheless, since the data are discrete, the most straightforward methodology is to discretize and approach the reconstruction problem as an algebraic problem (Section 2.5.6.4).
In digital image processing, space is represented by a multidimensional discrete lattice. (It is sometimes expedient to use cylindrical or spherical coordinates, but these also have to be appropriately discretized.) The 2D projections are sampled on a Cartesian grid , where d is the dimensionality of the grid (d = 2 for projections, d = 3 for the reconstructed object), is the size of the grid and a is the grid spacing (Fig. 2.5.6.4). In singleparticle reconstruction, the units of a are usually ångstroms and we also assume that the data were appropriately sampled, i.e., . Thus, the pixel size is less than or equal to the inverse of twice the maximum spatial frequency present in the data. Since the latter is not known in advance, a more practical rule is to select the pixel size at about one third of the expected resolution of the final structure, so that the adverse effects of interpolation are reduced.

Discretization in two dimensions (d = 2). The assumption is that the mass is located at the centre of the voxel. 
The input electron microscopy data (projections of the macromolecule) are discretized on a 2D Cartesian grid, but each projection has a particular orientation in polar coordinates. Except for a few cases (projection directions parallel to the axes of the coordinate system of the 3D structure), an interpolation is required to relate measured samples to the voxel (volume pixel) values on the 3D Cartesian grid of the reconstructed structure (Fig. 2.5.6.4). The step of backprojection can be visualized as a set of rays with base extended from projections and the ray values being added to the intersected voxels on the grid of the reconstructed structure (Figs. 2.5.6.3 and 2.5.6.4). One can select schemes that aim at approximation of the physical reality of the data collection, for example to weight the contributions by the areas of the voxels intersected by the ray or by the lengths of the lines that transverse the voxel (Huesman et al., 1977). In order to reduce the time of calculations, in electron microscopy one usually assumes that all the mass is located at the centre of the voxel and the additional accuracy is achieved by application of tri (or bi)linear interpolation. The exception is the algebraic reconstruction technique with blobs algorithm (Marabini et al., 1998), where the voxels are represented by smooth spherically symmetric volume elements [for example, the Keiser–Bessel function (2.5.6.44)].
In real space, both the projection and backprojection steps can be implemented in two different ways: as voxel driven or as ray driven (Laurette et al., 2000). If we consider a projection, in the voxeldriven approach the volume is scanned voxel by voxel. The nearest projection bin to the projection of each voxel is found, and the values in this bin and three neighbouring bins are increased by the corresponding voxel value multiplied by the weights calculated using bilinear interpolation. In the raydriven approach, the volume is scanned along the projection rays. The value of the projection bin is increased by the values in the volume calculated in equidistant steps along the rays using trilinear interpolation. Because voxel and raydriven methods apply interpolation to projections or to voxels, respectively, the interpolation artifacts will be different in each case. Therefore, when calculating reconstructions using iterative algorithms that alternate between projection and backprojection steps, it is important to maintain consistency; that is, to use the same method for both steps. In either case, the computational complexity of each method is of the order of K^{3}, although the voxeldriven approach is faster due to the smaller number of neighbouring points used in the interpolation.
In the reconstruction methods based on the direct Fourier inversion of the 3D ray transform, the interpolation is performed in Fourier space. Unfortunately, it is difficult to design an accurate and fast interpolation scheme for the discrete Fourier space. Bilinear interpolation introduces local errors and when applied in real space it results in attenuation of highfrequency information. When applied in Fourier space, bilinear interpolation results in errors evenly spread over the whole frequency range, thus resulting in potentially severe nonlocal errors in real space. In order to eliminate this error it would be tempting to use interpolation based on Shannon's sampling theorem [Shannon, 1949; reprinted in Proc. IEEE, (1998), 86, 447–457], which states that a properly sampled bandlimited signal can be fully recovered from its discrete samples. For the signal represented by K^{3} equispaced Fourier samples , the value of the Fourier transform at the arbitrary location (u, v, w) is given by (Crowther, Amos et al., 1970)where (Yuen & Fraser, 1979; Lanzavecchia & Bellon, 1994)In cryoEM, samples of are given at arbitrary 3D locations, as derived from Fourier transforms of 2D projection data (central section theorem) and one seeks to recover on the 3D Cartesian grid. Upon the inverse Fourier transform, it will yield the reconstructed object. The problem can be solved as an overdetermined system of linear equations (Crowther, DeRosier & Klug, 1970). Indeed, if we write and as 1D arrays and , respectively (the former has K^{2} × [number of projections] elements, while the latter has K^{3} elements), and we denote by W the appropriately dimensioned matrix of the interpolants, the leastsquares solution isThe above method is impractical because of the large size of the matrix W. In some cases, when due to symmetries the projection data are distributed approximately evenly (as in the case of icosahedral structures), the problem can be solved to a good degree of accuracy by performing the interpolation (2.5.6.17) independently along each of the three frequency axes (Crowther, DeRosier & Klug, 1970). Thus, in this case the solution to the problem of interpolation in Fourier space becomes a solution to the reconstruction problem.
For a more general singleparticle reconstruction application a moving Shannon window interpolation has been proposed (Lanzavecchia & Bellon, 1994, 1998). It is based on an attenuated version of the window function and in one dimension has the formwhere n is the window size and A is an integer that is even for K odd and vice versa. In multidimensional cases, a product of 's from (2.5.6.18) is used. In the case of interpolation between equispaced samples of , excellent results have been reported for n = 11 (Lanzavecchia et al., 1996). However, the application to the reconstruction problem, i.e., to resampling of the nonuniformly distributed Fourier data onto a 3D Cartesian grid, does not yield satisfactory results. Although general conditions under which interpolation using (2.5.6.18) can be done are known (Clark et al., 1985), they are not met in practice and the results are at best nonexact. In addition, the relatively large window size required (n = 11) results in impractical calculation times.
The algebraic methods have been derived based on the observation that when the projection equation (2.5.6.5) 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 and the voxels of the 3D object in a vector (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 K^{2}N and K^{3}, respectively, with K^{2} 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 are the interpolation weights. Their values are determined by the interpolation scheme used, but for the bi and trilinear interpolations . The algebraic version of (2.5.6.5) isMatrix P is rectangular and since in singleparticle reconstruction the number of projections exceeds the linear size of the object () the system of equations is overdetermined. It can be solved in a leastsquares sense:which yields a unique structure that corresponds to the minimum of . As in the case of direct inversion in Fourier space (2.5.6.17), the approach (2.5.6.20) is impractical because of the very large size of the projection matrix. Indeed, the size of P is K^{2}N × K^{3}, which for a modest number of projections N = 10 000 and image size K^{2} = 64^{2} = 4096 yields a matrix of size ~(4 × 10^{7}) × (2 × 10^{6})! 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) by using the singular value decomposition (SVD) of the matrix :where and contain eigenvectors of matrices and , respectively, and is a diagonal matrix containing the first r nonnegative square roots of the eigenvalues of (Golub & Van Loan, 1996). The solution is given bywhere is the inverse of the matrix containing the largest elements of 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).
In the general 3D case, a leastsquares 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 will contain only about 4N nonzero elements (out of the total ~K^{2}N). The main idea is that the matrix 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 ofby selecting the initial 3D structure (usually set to zero) and by iteratively updating its current approximation using the gradient of the objective function (2.5.6.23) ,Setting the relaxation parameter yields Richardson's algorithm (Gilbert, 1972). In singleparticle 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 singleparticle reconstruction (Frank et al., 1996) because it yields superior results over a wide range of experimental conditions (Penczek et al., 1992) 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 convergence by adjusting relaxation parameters: setting results in a steepest decent algorithm. Secondly, even faster convergence (in ~10 iterations) is achieved by solving (2.5.6.23) using the conjugategradient method, but this requires addition of a regularizing term to (2.5.6.23) in order to prevent excessive enhancement of noise. Such a term has the form , where matrix B is a discrete approximation of a Laplacian or higherorder 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) and by solving the problem for the structure g that in effect will be corrected for the CTF,Here we assumed that the projection data were grouped according to the defocus settings indexed by and we introduced a Lagrange multiplier whose value determines the smoothness of the solution g (Zhu et al., 1997). represents the spaceinvariant pointspread function of the microscope for the κth defocus group; thus has a blockToeplitz structure (Biemond et al., 1990) and . If we assume that there is no astigmatism (which means that the pointspread function is rotationally invariant), can be either applied to the projection of the structure or to the 3D structure itself. With these assumptions, (2.5.6.25) can be solved using the iterative scheme (Penczek et al., 1997)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)]. In addition, the product B^{T}B is the 3D Laplacian, which can be applied efficiently without actually creating the matrix B^{T}B. Finally, because of the large size of the matrix 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 realspace convolution. Algorithms (2.5.6.24)–(2.5.6.26) are implemented in the SPIDER package (Frank et al., 1996) and (2.5.6.24) in the SPARX package (Hohn et al., 2007).
The algebraic reconstruction technique (ART) predates SIRT; in the context of tomographic reconstructions it was proposed by Gordon and coworkers (Gordon et al., 1970) and later it was recognized as a version of Kaczmarz's method for iteratively solving (2.5.6.23) (Kaczmarz, 1993). We write (2.5.6.19) as a set of systems of equations, each relating single pixels f_{n}, in projections with the 3D structure,Note (2.5.6.27) and (2.5.6.19) are equivalent, as and . With this notation and a relaxation parameter , ART comprises the following steps:
Although the mathematical forms of update equations in SIRT (2.5.6.24) and in ART (2.5.6.28) 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) 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) in (2.5.6.28). 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; Herman & Meyer, 1993). Interestingly, a random order works almost equally well (Natterer & Wübbeling, 2001).
In singleparticle reconstruction, ART has been introduced in the form of `ART with blobs' (Marabini et al., 1998) and is available in the Xmipp package (Sorzano et al., 2004). 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, 1992; Matej & Lewitt, 1996). Introduction of blobs significantly reduces the number of iterations necessary to reach an acceptable solution (Marabini et al., 1998).
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; Youla & Webb, 1982; Sezan, 1992) introduced into EM by Carazo and coworkers (Carazo & Carrascosa, 1986, 1987; Carazo, 1992).
The method of filtered backprojection (FBP) is based on inversion formulae (2.5.6.11) (in two dimensions) or (2.5.6.14) (in three dimensions). It comprises the following steps: (i) a Fourier transform of each projection is computed; (ii) Fourier transforms of projections are multiplied by filters that account for a particular distribution of projections in Fourier space; (iii) the filtered projections are inversely Fourier transformed; (iv) realspace backprojection of processed projections yields the reconstruction. The method is particularly attractive due to the fact that the reconstruction calculated using simple realspace backprojection can be made efficient if the filter function is easy to compute.
In two dimensions with uniformly distributed projections the weighting function in Fourier space is the `ramp function' [(2.5.6.13)]. In two dimensions with nonuniformly distributed projections, when the analytical form of the distribution of projections is not known, an appropriate approximation to has to be found. A good choice is to select weights such that the backprojection integral becomes approximated by a Riemann sum (Penczek et al., 1996),For a given set of angles the weights are easily computed (Fig. 2.5.6.5).
In three dimensions, the weighting (2.5.6.29) is applicable in a singleaxis tilt datacollection geometry, where the 3D reconstruction can be calculated as a series of independent 2D reconstructions. In the general 3D case, the analogue of weighting (2.5.6.29) cannot be used, as the data are given in the form of 2D projections and it is not immediately apparent what fraction of the 3D Fourier volume is occupied by Fourier pixels in projections. However, the analogue of weighting (2.5.6.29) can be used in the inversion of 3D Radon transforms or in a direction inversion of 3D ray transforms that is based on an intermediate step of converting 2D projection data to 1D projection data, as described in Section 2.5.6.6.
In order to arrive at a workable solution, the weighting functions applicable to 2D projections are constructed based on an explicitly or implicitly formulated concept of the `local density' of projections. This concept was introduced by Bracewell (Bracewell & Riddle, 1967), who suggested for a 2D case of nonuniformly distributed projections a heuristic weighting function,The weighting function (2.5.6.30) can be easily extended to three dimensions; however, it has a major disadvantage that for a uniform distribution of projections it does not approximate well the weighting function (2.5.6.29), which we consider optimal.
Radermacher et al. (1986) proposed a derivation of a general weighting function using a deconvolution kernel calculated for a given (nonuniform) distribution of projections and, in modification of (2.5.6.14), a finite length of backprojection (Fig. 2.5.6.3). Such a `truncated' backprojection iswith projection directions andwhere D is the diameter of the object or the length of the backprojection l. By taking the Fourier transform of (2.5.6.31) and using the central section theorem (2.5.6.8), we obtain a 3D Fourier transform of the backprojected projection,and the 3D reconstruction is obtained by the inverse Fourier transform of the sum of contributions given by (2.5.6.33), and are variables in real and Fourier spaces, respectively, both extending in the direction of the projection direction .
In this analysis, the transfer function of the backprojection algorithm is obtained by setting , that is by finding the response of the algorithm to the input composed of delta functions in real space. This yields the inversion formula for a general weighted backprojection algorithm:with the weighting function given byThe general weighting function (2.5.6.36) is consistent with analytical solutions (2.5.6.11) and (2.5.6.14), as it can be shown that by assuming infinite support and continuous and uniform distribution of projection directions, in two dimensions one obtains (Radermacher, 2000).
The derivation of (2.5.6.36) is based on the analysis of continuous functions and its direct application to discrete data results in reconstruction artifacts; therefore, Radermacher (1992) proposed attenuating the sinc functions in (2.5.6.36) by exponent functions with decay depending on the diameter of the structure D, or simply replacing the sinc functions by exponent functions. This, however, reduces the concept of the weighting function corresponding to the deconvolution to the concept of the weighting function representing the `local density' of projections (2.5.6.30). The general weighted backprojection reconstruction is implemented in SPIDER (using exponentbased weighting functions) (Frank et al., 1996), Xmipp (Sorzano et al., 2004) and SPARX (Hohn et al., 2007).
Harauz & van Heel (1986) proposed basing the calculation of the density of projections, thus the weighting function, on the overlap of Fourier transforms of projections in Fourier space. Although the concept is general, it can be easily approached in two dimensions. If the diameter of the object is D, the width of a Fourier transform of a projection is 2/D (Fig. 2.5.6.5), as follows from the central section theorem (2.5.6.8). Harauz and van Heel postulated that the weighting should be inversely proportional to the sum of geometrical overlaps between a given central section and the remaining central sections. For a pair of projections il, this overlap iswhere T represents the triangle function (selected by the authors because it can be calculated efficiently). Also, owing to the Friedel symmetry of central sections, the angles in (2.5.6.37) are restricted such that . In this formulation, the overlap is limited to the maximum frequency,Thus, the overlap function becomesIn effect, the weighting function, called by the authors an `exact filter', is
The weighting function (2.5.6.40) easily extends to three dimensions; however, the calculation of the overlap between central sections in three dimensions (represented by slabs) is more elaborate (Harauz & van Heel, 1986). The method is conceptually simple and computationally efficient. However, (2.5.6.40) does not approximate the correct weighting well for a uniform distribution of projections [i.e., it should yield ]. This, as can be seen by integrating (2.5.6.39) over the whole angular range, is not the case. The exact filter backprojection reconstruction is implemented in the IMAGIC (van Heel et al., 1996), SPIDER (Frank et al., 1996) and SPARX (Hohn et al., 2007) packages.
The 3D reconstruction methods based on filtered backprojection are commonly used in singleparticle reconstruction. The reasons are: their versatility, ease of implementation, and – in comparison with iterative methods – good computational efficiency. Unlike in iterative methods, there are no parameters to adjust, although it has been noted that the results depend on the value of the diameter D of the structure in all three weighting functions [(2.5.6.30), (2.5.6.36) and (2.5.6.38)], so the performance of the reconstruction algorithm can be optimized for a particular datacollection geometry by changing the value of D (Paul et al., 2004). However, because computation of the weighting function involves calculation of pairwise distances between projections, the computational complexity is proportional to the square of the number of projections and for large data sets these methods become inefficient. It also has to be noted that the weighting functions (2.5.6.30), (2.5.6.36) and (2.5.6.40) remain approximations of the correct weighting function (2.5.6.29).
Direct Fourier methods are based on the central section theorem (2.5.6.8). A set of the 2D Fourier transforms of projections yields an approximation to on a nonuniform 3D grid, and a subsequent numerical 3D inverse Fourier transform gives an approximation to . If the 3D inverse Fourier transform could be realized by means of the 3D inverse fast Fourier transform (FFT), one would have a very fast reconstruction algorithm. Unfortunately, the preprocessing step yields on a nonuniform grid. In effect, the 3D inverse FFT is not applicable and an additional step of recovering samples of on a uniform grid from the available samples on a nonuniform grid is necessary.
One possibility is to resample the nonuniformly sampled version of onto a 3D Cartesian grid by some form of interpolation. For example, Grigorieff used a modified trilinear interpolation scheme in the FREALIGN package (Grigorieff, 1998). Simple interpolation methods have been found to give inaccurate results, although more sophisticated interpolation schemes can go a long way to improve the accuracy (Lanzavecchia et al., 1993). Unfortunately, the recommended window size makes them impractical for most applications.
The most accurate Fourier reconstruction methods are those that employ nonuniform Fourier transforms, particularly the 3D gridding method (O'Sullivan, 1985; Schomberg & Timmer, 1995). The griddingbased direct Fourier reconstruction algorithm (GDFR) (Penczek et al., 2004) was developed specifically for singleparticle reconstruction. It comprises three steps:
The method involves a number of parameters. First, we need to decide the oversampling factor for the padding with zeros of input projections before FFTs are computed. In GDFR, the oversampling is set to two, although smaller factors can also yield good results. Second, we need a window function whose support in Fourier space is `small'. In order to assure good computational efficiency of the convolution step (2.5.6.41), in GDFR this support was set to six Fourier voxels. In addition, in order to prevent division by zeros in (2.5.6.43), the weighting function must be positive within the support of the reconstructed object. A recommendable window function is the separable, bellshaped Kaiser–Bessel window (O'Sullivan, 1985; Jackson et al., 1991; Schomberg & Timmer, 1995) (Fig. 2.5.6.6a):where is the zeroorder modified Bessel function and α = 1.25 is a parameter. The weighting function associated with isFinally, the gridding weights c are chosen such that in discrete implementation of (2.5.6.41) and (2.5.6.42) we obtain a Riemann sum approximating the respective integral (Schomberg, 2002). The Voronoi diagram (Okabe et al., 2000) of the sampling points provides a good partitioning of sampling space, so the gridding weights c are computed by constructing a Voronoi diagram of the grid points and by choosing the weights as the volumes of the Voronoi cells (Fig. 2.5.6.7a). In cryoEM, the number of projections, thus the number of sampling points in Fourier space, is extremely large. Thus, although the calculation of the gridding weights via the 3D Voronoi diagram of the nonuniformly spaced grid points for would lead to an accurate direct Fourier method, the method would be very slow and would require excessive computer memory. To circumvent this problem, in GDFR (Penczek et al., 2004) the 2D reverse gridding method is used to compute the Fourier transform of each projection on a 2D polar grid. In this way, is obtained on a 3D spherical grid, where the grid points are located both on centred spheres and on straight lines through the origin. Accordingly, it becomes possible to partition the sampled region into suitable sampling cells via the computation of a 2D Voronoi diagram on a unit sphere, rather than a 3D Voronoi diagram in Euclidean space (Fig. 2.5.6.7b). This significantly reduces the memory requirements and improves the computational efficiency of the method, particularly when a fast algorithm for computing the Voronoi diagram on the unit sphere is employed (Renka, 1997).
The gridding weights in GDFR correspond to the weighting functions in filtered backprojection algorithms. Moreover, setting the gridding weights to be equal to angular regions on the unit sphere directly corresponds to the equating of weighting function (2.5.6.29) in 2D filtered backprojection to the length of an arc on a unitary circle. Thus, both methods yield nearly optimum projection density functions.
The reversed gridding method is obtained by reversing the sequence of steps (1)–(3) of the gridding method:
Note the reverse gridding does not require explicit gridding weights, as in the third step they are constant (Penczek et al., 2004).
In GDFR, the third step of the reversed gridding results in a set of 1D Fourier central lines calculated using a constant angular step. Clearly, upon inverse Fourier transform they amount to a Radon transform of the projection and, upon repeating the process for all available projections, they yield a Radon transform of a 3D object (albeit nonuniformly sampled with respect to Eulerian angles). Thus GDFR, in addition to being a method of inverting a 3D ray transform, is also a highly accurate method of inverting a 3D Radon transform. GDFR is implemented in the SPIDER package (Frank et al., 1996).
The results of a comparison of selected reconstruction algorithms are shown in Fig. 2.5.6.8. The tests were performed using simulated data with the projections of a phantom 3D structure calculated using the inverse gridding method (Penczek et al., 2004). The GDFR yields a virtually perfect reconstruction that agrees with the phantom over the whole frequency range (with the exception of the highest frequencies, but these cannot be reproduced due to geometrical limitations). GD3D is also a griddingbased reconstruction algorithm, in which the gridding weights are calculated directly from contributions of the weighting function to the 3D Fourier pixels and are applied to the voxels in the 3D Fourier volume (Jackson et al., 1991). This yields an approximation to the `local density' of contributing projections. The simplified gridding weights in GD3D result in deterioration of the reconstruction in the low and intermediatefrequency ranges. The general weighted backprojection with exponentbased weighting function (WBP1) performs quite well, although reproduction in the lowfrequency range is inferior. Similar artifacts are present in the reconstruction using the exact filter weighted backprojection (WBP2), which, in addition, performs disappointingly at higher frequencies. This relatively poor performance is attributed to the nonoptimal weighting schemes used in both methods. It is also interesting to note that the backprojection step is identical in both algorithms; they only differ by the weighting function. The significant difference in their performance attests to the importance of good weighting schemes for highquality 3D reconstructions from nonuniformly distributed projections. SIRT yields a reconstruction that in the low and intermediatefrequency ranges matches in quality the reconstruction obtained with GDFR. However, there is a significant loss of quality at high frequencies. This is due to the inferior linear interpolation scheme used in SIRT (this impediment is shared with WBP1 and WBP2). [The dip of the SIRT fidelity curve at a spatial frequency of 0.42 in Fig. 2.5.6.8(a) is due to an inconsistency between the interpolation methods used internally in SIRT to generate intermediate projections of the object to be reconstructed and the Fourierspacebased method used to generate the test data.] In terms of computational efficiency, griddingbased algorithms outperform weighted backprojection algorithms by a small factor while SIRT is approximately ten times slower (depending on the number of iterations used).
Many objects imaged by EM have symmetries; the two types that are often met are helical (phage tails, Factin, microtubules, myosin thick filaments, and bacterial pili and flagella) and pointgroup symmetries (many macromolecular assemblies, virus capsids). If the object has helical symmetry, it is convenient to use cylindrical coordinates and a dedicated reconstruction algorithm; particularly for the initial analysis of the data that has to be done in Fourier space. Diffraction from such structures with c periodicity and scattering density is defined by the Fourier–Bessel transform:The inverse transform has the form so that and are the mutual Bessel transforms and
Owing to helical symmetry, (2.5.6.48) and (2.5.6.49) contain only those of the Bessel functions that satisfy the selection rule (Cochran et al., 1952) where N, q and p are the helix symmetry parameters, . Each layer l is practically determined by the single function with the lowest n; the contributions of other functions are neglected. Thus, the Fourier transform of one projection of a helical structure, with an account of symmetry and phases, gives the 3D transform (2.5.6.49). However, biological specimens tend to be flexible and disordered, and exact helical symmetry is rarely observed. A possible approach to dealing with flexibility is to computationally straighten filaments (Egelman, 1986), but this has the potential for introducing artifacts. Another difficulty with helical analysis is that the indexing of a pattern can be ambiguous and the wrong symmetry can be chosen (Egelman & Stasiak, 1988). Further complications exist when the filament does not have a precisely defined helical symmetry, such as Factin, which has a random variable twist (Egelman et al., 1982). To address these problems, Egelman developed a realspace refinement method for the reconstruction of helical filaments that is capable of determining the helical symmetry of an unknown structure (Egelman, 2000). In this approach, the Fourier–Bessel reconstruction algorithm is replaced by a general reconstruction algorithm discussed in Section 2.5.6.6 with the realspace projections supplied as input. Thus, the symmetry is not enforced within the reconstruction algorithm; instead, it is determined and imposed subsequently by realspace averaging.
The presence of pointgroup symmetry in the structure means that any projection yields as many symmetryrelated (but differently oriented) copies of itself as the number of symmetry operations in the group. For example, for three dimensions, each projection enters the reconstruction process at six different projection directions, while for icosahedral symmetry (I) the number is sixty! This high multiplicity makes it possible to write a reconstruction program that explicitly takes into account the symmetry and can perform the task much faster than generic algorithms. Such programs are usually part of dedicated software packages that were specifically designed for the determination of icosahedral structures (Fuller et al., 1996; Lawton & Prasad, 1996; Liang et al., 2002). These icosahedral reconstruction programs are heavily indebted to a set of FORTRAN routines written by R. A. Crowther at the Medical Research Council Laboratory of Molecular Biology (MRC LMB) and perform the 3D reconstruction using the Fourier–Bessel transform strategy outlined above [(2.5.6.46)–(2.5.6.49)] (Crowther, Amos et al., 1970; Crowther & Amos, 1972).
The general reconstruction methods (algebraic, filtered backprojection, direct Fourier inversion) easily accommodate symmetries in the data. In major singleparticle reconstruction packages, reconstruction programs are implemented such that the pointgroup symmetry is a parameter of the program. The symmetry operation is internally taken into account during calculation of the weighting function and it is also applied to the set of Eulerian angles assigned to each projection so the multiple copies are implicitly created and processed. This approach results in an extended time of calculations, but it is entirely general. In direct Fourier inversion algorithms the numerical inaccuracy of the symmetrization performed in Fourier space will result in nonsymmetric artifacts in real space. Thus, in SPIDER (Frank et al., 1996) an additional realspace symmetrization is performed after the reconstruction is completed. Finally, it has to be noted that although it might be tempting to calculate a 3D reconstruction without enforcement of symmetry and to symmetrize the resulting structure, this approach is incorrect. This fact can be seen from the way weighting functions are constructed: weighting functions calculated with symmetries taken into account are not equal to the weighting function calculated for a unique set of projections and subsequently symmetrized.
References
Biemond, J., Lagendijk, R. L. & Mersereau, R. M. (1990). Iterative methods for image deblurring. Proc. IEEE, 78, 856–883.Bracewell, R. N. (1956). Strip integration in radio astronomy. Austr. J. Phys. 9, 198–217.
Bracewell, R. N. & Riddle, A. C. (1967) Inversion of fanbeam scans in radio astronomy. Astrophys. J. 150, 427–434.
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 threedimensional reconstructions of crystalline specimens by the method of convex projections. J. Microsc. 145, 159–177.
Clark, J. J., Palmer, M. R. & Lawrence, P. D. (1985). A transformation method for the reconstruction of functions from nonuniformly spaced samples. IEEE Trans. Acoust. Speech Signal Process. 33, 1151–1165.
Cochran, W., Crick, F. H. C. & Vand, V. (1952). The structure of synthetic polypeptides. 1. The transform of atoms on a helix. Acta Cryst. 5, 581–586.
Crowther, R. A. & Amos, L. A. (1972). Threedimensional image reconstructions of some small spherical viruses. Cold Spring Harbor Symp. Quant. Biol. 36, 489–494.
Crowther, R. A., Amos, L. A., Finch, J. T., DeRosier, D. J. & Klug, A. (1970). Three dimensional reconstruction of spherical viruses by Fourier synthesis from electron micrographs. Nature (London), 226, 421–425.
Crowther, R. A., DeRosier, D. J. & Klug, A. (1970). The reconstruction of a threedimensional structure from projections and its application to electron microscopy. Proc. R. Soc. London Ser. A, 317, 319–340.
DeRosier, D. J. & Klug, A. (1968). Reconstruction of three dimensional structures from electron micrographs. Nature (London), 217, 130–134.
Egelman, E. (1986). An algorithm for straightening images of curved filamentous structures. Ultramicroscopy, 19, 367–374.
Egelman, E. H. (2000). A robust algorithm for the reconstruction of helical filaments using singleparticle methods. Ultramicroscopy, 85, 225–234.
Egelman, E. H., Francis, N. & DeRosier, D. J. (1982). Factin is a helix with a random variable twist. Nature (London), 298, 131–135.
Egelman, E. H. & Stasiak, A. (1988). Structure of helical RecA–DNA complexes. II. Local conformational changes visualized in bundles of RecA–ATP gamma S filaments. J. Mol. Biol. 200, 329–349.
Frank, J. (2006). ThreeDimensional Electron Microscopy of Macromolecular Assemblies. New York: Oxford University Press.
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.
Fuller, S. D., Butcher, S. J., Cheng, R. H. & Baker, T. S. (1996). Threedimensional reconstruction of icosahedral particles – the uncommon line. J. Struct. Biol. 116, 48–55.
Gilbert, P. F. C. (1972). Iterative methods for the threedimensional 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 threedimensional electron microscopy and Xray photography. J. Theor. Biol. 29, 471–481.
Grigorieff, N. (1998). Threedimensional structure of bovine NADH: ubiquinone oxidoreductase (complex I) at 22 Å in ice. J. Mol. Biol. 277, 1033–1046.
Hamaker, C. & Solmon, D. C. (1978). Angles between null spaces of Xrays. J. Math. Anal. Appl. 62, 1–23.
Harauz, G. & van Heel, M. (1986). Exact filters for general geometry threedimensional reconstruction. Optik, 73, 146–156.
Heel, M. van, Harauz, G. & Orlova, E. V. (1996). A new generation of the IMAGIC image processing system. J. Struct. Biol. 116, 17–24.
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 cryoEM image processing. J. Struct. Biol. 157, 47–55.
Holmes, K. C., Angert, I., Kull, F. J., Jahn, W. & Schroder, R. R. (2003). Electron cryomicroscopy shows how strong binding of myosin to actin releases nucleotide. Nature (London), 425, 423–427.
Huesman, R. H., Gullberg, G. T., Greenberg, W. L. & Budinger, T. F. (1977). RECLBL library users' manual – Donner algorithms for reconstruction tomography. University of California, Berkeley, USA.
Jackson, J. I., Meyer, C. H., Nishimura, D. G. & Macovski, A. (1991). Selection of a convolution function for Fourier inversion using gridding. IEEE Trans. Med. Imaging, 10, 473–478.
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.
Lanzavecchia, S. & Bellon, P. L. (1994). A moving window Shannon reconstruction for image interpolation. J. Visual Comm. Image Repres. 3, 255–264.
Lanzavecchia, S. & Bellon, P. L. (1998). Fast computation of 3D Radon transform via a direct Fourier method. Bioinformatics, 14, 212–216.
Lanzavecchia, S., Bellon, P. L. & Scatturin, V. (1993). SPARK, a kernel of software programs for spatial reconstruction in electron microscopy. J. Microsc. 171, 255–266.
Lanzavecchia, S., Tosoni, L. & Bellon, P. L. (1996). Fast sinogram computation and the sinogrambased alignment of images. Comp. Appl. Bio. Sci. 12, 531–537.
Laurette, I., Zeng, G. L., Welch, A., Christian, P. E. & Gullberg, G. T. (2000). A threedimensional raydriven attenuation, scatter and geometric response correction technique for SPECT in inhomogeneous media. Phys. Med. Biol. 45, 3459–3480.
Lawton, J. A. & Prasad, B. V. V. (1996). Automated software package for icosahedral virus reconstruction. J. Struct. Biol. 116, 209–215.
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.
Liang, Y. Y., Ke, E. Y. & Zhou, Z. H. (2002). IMIRS: a highresolution 3D reconstruction package integrated with a relational image database. J. Struct. Biol. 137, 292–304.
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 3D 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.
O'Sullivan, J. D. (1985). A fast sinc function gridding algorithm for Fourier inversion in computer tomography. IEEE Trans. Med. Imaging, 4, 200–207.
Okabe, A., Boots, B., Sugihara, K. & Chiu, S. N. (2000). Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. New York: John Wiley & Sons.
Orlov, S. S. (1976). Theory of threedimensional reconstruction 1. Conditions for a complete set of projections. Sov. Phys. Crystallogr. 20, 312–314.
Paul, D., Patwardhan, A., Squire, J. M. & Morris, E. P. (2004) Single particle analysis of filamentous and highly elongated macromolecular assemblies. J. Struct. Biol. 148, 236–250.
Penczek, P., Radermacher, M. & Frank, J. (1992). Threedimensional reconstruction of single particles embedded in ice. Ultramicroscopy, 40, 33–53.
Penczek, P. A. (2002). Threedimensional spectral signaltonoise ratio for a class of reconstruction algorithms. J. Struct. Biol. 138, 34–46.
Penczek, P. A., Renka, R. & Schomberg, H. (2004). Griddingbased direct Fourier inversion of the threedimensional ray transform. J. Opt. Soc. Am. A, 21, 499–509.
Penczek, P. A., Zhu, J. & Frank, J. (1996). A commonlines based method for determining orientations for N > 3 particle projections simultaneously. Ultramicroscopy, 63, 205–218.
Penczek, P. A., Zhu, J., Schröder, R. & Frank, J. (1997). Threedimensional reconstruction with contrast transfer function compensation from defocus series. Scan. Microsc. Suppl. 11, 1–10.
Radermacher, M. (1992). Weighted backprojection methods. In Electron Tomography, edited by J. Frank, pp. 91–115. New York: Plenum.
Radermacher, M. (1994). Threedimensional reconstruction from random projections: orientational alignment via Radon transforms. Ultramicroscopy, 53, 121–136.
Radermacher, M. (2000). Threedimensional reconstruction of single particles in electron microscopy. In Image Analysis: Methods and Applications, edited by D.P. Häder, pp. 295–328. Boca Raton: CRC Press.
Radermacher, M., Wagenknecht, T., Verschoor, A. & Frank, J. (1986). A new 3D reconstruction scheme applied to the 50S ribosomal subunit of E. coli. J. Microsc. 141, RP1–RP2.
Renka, R. J. (1997). Algorithm 772. STRIPACK: Delaunay triangulation and Voronoi diagram on the surface of a sphere. ACM Trans. Math. Software, 23, 416–434.
Schomberg, H. (2002). Notes on direct and griddingbased Fourier inversion methods. In Proceedings of the IEEE International Symposium on Biomedical Imaging, Washington, DC, edited by M. Unser & Z. P. Liang, pp. 645–648.
Schomberg, H. & Timmer, J. (1995). The gridding method for image reconstruction by Fourier transformation. IEEE Trans. Med. Imaging, 14, 596–607.
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.
Shannon, C. E. (1949). Communication in the presence of noise. Proc. Inst. Radio Eng. 37, 10–21.
Sorzano, C. O. S., Marabini, R., VelazquezMuriel, J., BilbaoCastro, J. R., Scheres, S. H. W., Carazo, J. M. & PascualMontano, A. (2004). XMIPP: a new generation of an opensource image processing package for electron microscopy. J. Struct. Biol. 148, 194–204.
Vainshtein, B. K. (1971). Finding the structure of objects from projections. Sov. Phys. Crystallogr. 15, 781–787.
Youla, D. C. & Webb, H. (1982). Image restoration by the method of convex projections. 1. Theory. IEEE Trans. Med. Imaging, 1, 81–94.
Yuen, C. K. & Fraser, D. (1979). Digital Spectral Analysis. Adelaide: CSIRO Pitman.
Zhu, J., Penczek, P. A., Schröder, R. & Frank, J. (1997). Threedimensional reconstruction with contrast transfer function correction from energyfiltered cryoelectron micrographs: procedure and application to the 70S Escherichia coli ribosome. J. Struct. Biol. 118, 197–219.