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

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

Section 2.5.6. Three-dimensional reconstruction5

B. K. Vainshteinc and P. A. Penczekg

2.5.6. Three-dimensional reconstruction5

| top | pdf |

2.5.6.1. The object and its projection

| top | pdf |

In electron microscopy (EM) and single-particle reconstruction, three-dimensional (3D) reconstruction methods are used for studying biological structures; that is, symmetric or asymmetric associations of biomacromolecules (muscles, spherical and rod-like viruses, bacteriophages, individual proteins and ribosomes) (Frank, 2006[link]). The electron microscope is used to obtain parallel-beam two-dimensional (2D) projections [\varphi _2 ( {\bf x},\boldtau )] (ray transform) of frozen hydrated 3D macromolecules [\varphi _3 ( {\bf r})] suspended in random orientations (Fig. 2.5.6.1[link]). The function [\varphi _2( {\bf x}_{\boldtau })] is the 2D projection of the 3D molecular electron distribution [\varphi _3 ( {\bf r})]. One can also consider one-dimensional (1D) projections [\varphi _1 ( s,\boldtau )] 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[link]).

[Figure 2.5.6.1]

Figure 2.5.6.1 | top | pdf |

A three-dimensional object [\varphi_{3}] and its two-dimensional projection [\varphi_{2}]. The electron beam penetrates the specimen in the direction of the z axis.

Within the linear weak-phase-object approximation of the image-formation process in the microscope [see equation (2.5.2.42)[link] in Section 2.5.2[link] of this chapter], 2D projections represent line integrals of the potential of the particle under examination convoluted with the point-spread function of the microscope, s, so, using (2.5.2.43)[link],[I(xy) = 1 + 2\sigma s(xy) * \varphi (xy) = 1 + 2\sigma s(xy) * \textstyle \int \varphi _3 ( {\bf r})\,{\rm d}z.\eqno(2.5.6.1)]Since[\varphi_{2} ({\bf x}_{\boldtau}) = {\textstyle\int} \varphi_{3} ({\bf r})\ {\rm d}\tau, \quad {\boldtau} \perp {\bf x,} \eqno(2.5.6.2)]we have[I(xy) = 1 + 2\sigma s( {\bf x}_\boldtau ) * \varphi _2 ( {\bf x}_\boldtau),\quad \boldtau \perp {\bf x}. \eqno(2.5.6.3)]If we omit constant terms, we obtain[I(xy) = s( {\bf x}_\boldtau) * \varphi _2 ( {\bf x}_\boldtau) = s( {\bf x}_\boldtau ) * \textstyle \int \varphi _3 ( {\bf r})\, {\rm d}\tau, \quad \boldtau \perp {\bf x}. \eqno(2.5.6.4)]

In this section, we will assume that all images were collected using the same defocus setting, so the point-spread function s is constant and does not depend on the projection direction [\boldtau]. Thus, we will concern ourselves with the inversion of the projection problem[\varphi _2 ( {\bf x}_\boldtau ) = \textstyle \int \varphi _3 ( {\bf r})\,{\rm d}\tau,\quad \boldtau \perp {\bf x}. \eqno(2.5.6.5)]The projection direction is defined by a unit vector [\boldtau ( \theta ,\psi )] and it is formed on the plane x perpendicular to [\boldtau]. The set of various projections [\varphi _2 ({\bf x}_{\boldtau_i }) = \varphi _{2_i }( {\bf x}_i )] may be assigned by a discrete or continuous set of points on a unit sphere [| \tau | = 1] (Fig. 2.5.6.2[link]).

[Figure 2.5.6.2]

Figure 2.5.6.2 | top | pdf |

The projection sphere and projection [\varphi _2 ( {{\bf x}_\boldtau } )] of [\varphi_{3}({\bf r})] along τ onto the plane [{\boldtau} \perp {\bf x}]. The case [{\boldtau} \perp Z] 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 [\varphi _2 ] of a 3D object [\varphi _3 ] 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[link]; DeRosier & Klug, 1968[link]; Crowther, Amos et al., 1970[link]). In Cartesian coordinates, a 3D Fourier transform is[ \eqalignno{ {\scr F}[\varphi_{3} ({\bf r})] &= \Phi_{3} (u,v,w) &\cr &= {\textstyle\int\!\!\int\!\!\int} \varphi_{3} (x,y,z) \exp \{2 \pi i (ux + vy + wz)\} \ \hbox{d}x \ \hbox{d}y \ \hbox{d}z. &\cr &&(2.5.6.6)\cr}]The transform of projection [\varphi _2 (x,y)] along z is[ \eqalignno{ {\scr F} [\varphi_{2} (x,y)] &= \Phi_{3} (u,v,0) &\cr &= {\textstyle\int\!\!\int\!\!\int} \varphi_{3} (x,y,z) \exp \{2 \pi i (ux + vy + 0z)\} \ \hbox{d}x \ \hbox{d}y \ \hbox{d}z &\cr &= {\textstyle\int\!\!\int\!\!\int} \varphi_{3} (x,y,z) \ \hbox{d}z \exp \{2 \pi i (ux + vy)\} \ \hbox{d}x \ \hbox{d}y &\cr &= {\textstyle\int\!\!\int} \varphi_{2} (x,y) \exp \{2 \pi i (ux + vy)\} \ \hbox{d}x \ \hbox{d}y &\cr &= \Phi_{2} (u,v). &(2.5.6.7)\cr}]In the general case of projecting along the vector [\boldtau], the central section theorem is[ {\scr F}[\varphi_{2} ({\bf x}_{\boldtau})] = \Phi_{3} ({\bf u}_{\boldtau}),\quad {\bf u}_\boldtau \perp \boldtau. \eqno(2.5.6.8)]From this theorem it follows that the inversion of the 3D ray transform is possible if there is a continuous set of projections [\varphi _\boldtau ] corresponding to the motion of the vector [\boldtau ( \theta ,\psi )] over any continuous line connecting the opposite points on the unit sphere (Fig. 2.5.6.2[link]) (Orlov's condition: Orlov, 1976[link]). This result is evidenced by the fact that in this case the cross sections [{\scr F}[ \varphi _2 ( {\bf x}_{\boldtau } )]] that are perpendicular to [\boldtau ] in Fourier space continuously cover the whole Fourier space, i.e., they yield [{\scr F}[ \varphi _3( {\bf r} )]] and thereby determine [\varphi _3 ( {\bf r}) = {\scr F}^{ - 1} [ \Phi _3 ( {\bf u})]].

In single-particle reconstruction, imaged objects are randomly and nonuniformly oriented on the substrate at different angles (Frank, 2006[link]) 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 signal-to-noise ratio (Penczek, 2002[link]).

2.5.6.2. 3D reconstruction in the general case

| top | pdf |

In the general case of the 3D reconstruction of [\varphi _3 ( {\bf r})] from projections [\varphi _2 ( {\bf x}_\boldtau)], the projection vector [\boldtau] occupies arbitrary positions on the projection sphere (Fig. 2.5.6.2[link]). First, let us consider the case of a 2D function [\rho _2 ( {\bf x})] and its ray transform [\varphi _1 ( x,\psi )]. We introduce an operation of backprojection b, which is stretching along [\boldtau _\psi ] each 1D projection [\varphi _1 ( x_\psi )] (Fig. 2.5.6.3[link]). When the result is integrated over the full angular range of projections [\varphi _1 ( x,\psi )], we obtain the projection synthesis defined as[b(x,y ) = \textstyle \int \limits_0^\pi {\varphi _1 ( {x\cos \psi + y\sin \psi ,\psi } )} \,{\rm d}\psi. \eqno(2.5.6.9)]However, the backprojection operator is not the inverse of a 2D ray transform, as the resulting image b is blurred by the point-spread function [( x^2 + y^2 )^{-1/2}] (Vainshtein, 1971[link]):[b(x,y) = \rho ( x,y ) * ( {x^2 + y^2 } )^{-1/2}.\eqno(2.5.6.10)]By noting that the Fourier transform of [( x^2 + y^2 )^{-1/2}] is [( u^2\ +] [ v^2 )^{ - 1/2}] and by using the convolution theorem [{\scr F}[ f * g] =] [ {\scr F}[ f ]{\scr F}[ g ]], we obtain the `backprojection-filtering' inversion formula:[\eqalignno{\rho ( x,y ) &= b( x,y ) * ( x^2 + y^2 )^{1/2} = {\scr F}^{ - 1}\big[ | {\bf u} |{\scr F}[ b ] \big] &\cr& ={\rm Filtration}_{| {\bf u} |} [ {\rm Backprojection}( \varphi _1 ) ].&(2.5.6.11)}]The more commonly used `filtered-backprojection' inversion is based on the 2D version of the central section theorem (2.5.6.8)[link]:[{\scr F}[ \varphi _1 ( x_\psi ) ] = \Upsilon _2 ( {\bf u}_\psi ) = \Upsilon _2 ( R,\Psi ),\eqno(2.5.6.12)]where [{\scr F}[ \rho _2 ] = \Upsilon _2 ]. With this in mind, [\rho _2 ( {\bf x} )] can be related to its ray transform by evaluating the Fourier transform of [\rho _2 ] in polar coordinates:[\eqalignno{ \rho _2 ( {\bf x} ) &= \textstyle \int \Upsilon _2 ( {\bf u} )\exp ( - 2\pi i{\bf ux} )\, {\rm d}{\bf u} &\cr & = \textstyle \int \limits_0^\pi \textstyle \int \limits_{ - \infty }^\infty \Upsilon _2 ( R,\Psi )\exp ( - 2\pi i{\bf x}\boldtau )| R |\,{\rm d}R\,{\rm d}\Psi&\cr & = \textstyle \int \limits_0^\pi \textstyle \int \limits_{ - \infty }^\infty {\scr F}[ \varphi _1 ( x_\psi ) ]\exp ( - 2\pi i{\bf x}\boldtau )| R |\,{\rm d}R\,{\rm d}\Psi &\cr & =\textstyle \int \limits_0^\pi {\scr F}^{ - 1}\big [ | R |{\scr F}[ \varphi _1 ( x_\psi ) ] \big]\,{\rm d}\Psi &\cr & = {\rm Backprojection}[ {\rm Filtration}_{| R |} ( \varphi _1 )]. &(2.5.6.13)} ]In three dimensions, the backprojection stretches each 2D projection [\varphi _{2_i} [ {\bf x},{\boldtau }( \theta ,\psi )_i ]] along the projection direction [\boldtau( \theta ,\psi )_i ]. A 3D synthesis is the integral over the hemisphere (Fig. 2.5.6.2[link])[b( {\bf r} ) = \textstyle \int \limits_\omega \varphi _2 ( {\bf x},\omega _\boldtau )\,{\rm d}\omega _\boldtau = \varphi _3 ( {\bf r} ) * ( x^2 + y^2 + z^2 )^{ - 1}. \eqno(2.5.6.14)]Thus, in three dimensions the image b obtained using the backprojection operator is blurred by the point-spread function [1/( x^2 + y^2 + z^2 )] (Vainshtein, 1971[link]). It is possible to derive inversion formulae analogous to (2.5.6.11)[link] and (2.5.6.13)[link].

[Figure 2.5.6.3]

Figure 2.5.6.3 | top | pdf |

(a) Formation of a backprojection function; (b) projection synthesis (2.5.6.9)[link] 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)[link]. 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)[link] and (2.5.6.14)[link], a simple backprojection step results in reconstruction that corresponds to a convolution of the original function with a point-spread 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 point-spread 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[link]). Similarly, the central section theorem gives rise to direct Fourier inversion algorithms (Section 2.5.6.6[link]). 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[link]).

2.5.6.3. Discretization and interpolation

| top | pdf |

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 [\varphi _2 ( {\bf x} )] are sampled on a Cartesian grid [\{ {\bf k}a:{\bf k} \in {\bf Z}^d ,(-{\bf K}/2) \leq{\bf k}\,\le\,({\bf K}/2)\}], where d is the dimensionality of the grid (d = 2 for projections, d = 3 for the reconstructed object), [{\bf K} \in {\bf Z}_ + ^d ] is the size of the grid and a is the grid spacing (Fig. 2.5.6.4[link]). In single-particle reconstruction, the units of a are usually ångstroms and we also assume that the data were appropriately sampled, i.e., [a \leq (1/2u_{\rm max})]. 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.

[Figure 2.5.6.4]

Figure 2.5.6.4 | top | pdf |

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[link]). The step of backprojection can be visualized as a set of rays with base [a^{d - 1} ] 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[link] and 2.5.6.4[link]). 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[link]). 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[link]), where the voxels are represented by smooth spherically symmetric volume elements [for example, the Keiser–Bessel function (2.5.6.44)[link]].

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[link]). If we consider a projection, in the voxel-driven 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 ray-driven 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 ray-driven 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 K3, although the voxel-driven 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 high-frequency 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[link]; reprinted in Proc. IEEE, (1998), 86, 447–457], which states that a properly sampled band-limited signal can be fully recovered from its discrete samples. For the signal represented by K3 equispaced Fourier samples [\Phi _{3hkl} ], the value of the Fourier transform [\Phi _3 ] at the arbitrary location (u, v, w) is given by (Crowther, Amos et al., 1970[link])[\Phi _3 ( {u,v,w} ) = \textstyle\sum\limits_{h = 0}^{K - 1} {\textstyle\sum\limits_{k = 0}^{K - 1} {\textstyle\sum\limits_{l = 0}^{K - 1} {\Phi _{3hkl} w_h ( u )} } } w_k ( v )w_l ( w ), \eqno(2.5.6.15)]where (Yuen & Fraser, 1979[link]; Lanzavecchia & Bellon, 1994[link])[\let\normalbaselines\relax\openup 5pt w_k ( u ) =\left \{ \matrix{{\displaystyle\sin\{K\pi [u-(k/K)]\}\over \displaystyle\sin\{\pi [u-(k/K)]\}},\hfill & K\,\,{\rm odd}\semi\hfill \cr {\displaystyle\sin\{K\pi [u-(k/K)]\}\over \displaystyle\tan\{\pi [u-(k/K)]\}},\hfill &K\,\,{\rm even}.\hfill}\right.\eqno(2.5.6.16)]In cryo-EM, samples of [\Phi _3 ] are given at arbitrary 3D locations, as derived from Fourier transforms of 2D projection data (central section theorem) and one seeks to recover [\Phi _{3hkl} ] 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[link]). Indeed, if we write [\Phi _3 ] and [\Phi _{3hkl} ] as 1D arrays [{\boldPhi }_3 ] and [{\boldPhi }_{3( {hkl} )} ], respectively (the former has K2 × [number of projections] elements, while the latter has K3 elements), and we denote by W the appropriately dimensioned matrix of the interpolants, the least-squares solution is[{\boldPhi }_{3( {hkl} )} = ( {{\bf{W}}^T {\bf{W}}} )^{ - 1} {\bf{W}}^T {\boldPhi }_3.\eqno(2.5.6.17)]The 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)[link] independently along each of the three frequency axes (Crowther, DeRosier & Klug, 1970[link]). 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 single-particle reconstruction application a moving Shannon window interpolation has been proposed (Lanzavecchia & Bellon, 1994[link], 1998[link]). It is based on an attenuated version of the window function and in one dimension has the form[\eqalignno{\Phi _1 ( u ) &= \displaystyle\sum\limits_{k = m}^{m + n - 1} \Phi _{1k} {\sin \{\pi [u-(k/K)]\}\over \sin\{(\pi/2)[u-(k/K)]\}}&\cr &\quad\times \left(\cos\{(\pi/2)[u-(k/K)]\}\right)^A,&(2.5.6.18)}]where n is the window size [( {n \ll K} )] and A is an integer that is even for K odd and vice versa. In multidimensional cases, a product of [\Phi]'s from (2.5.6.18)[link] is used. In the case of interpolation between equispaced samples of [\Phi _{2hk} ], excellent results have been reported for n = 11 (Lanzavecchia et al., 1996[link]). 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)[link] can be done are known (Clark et al., 1985[link]), 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.

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]).

2.5.6.5. Filtered backprojection

| top | pdf |

The method of filtered backprojection (FBP) is based on inversion formulae (2.5.6.11)[link] (in two dimensions) or (2.5.6.14)[link] (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) real-space backprojection of processed projections yields the reconstruction. The method is particularly attractive due to the fact that the reconstruction calculated using simple real-space backprojection can be made efficient if the filter function is easy to compute.

In two dimensions with uniformly distributed projections the weighting function [c(| R |,\Psi )] in Fourier space is the `ramp function' [| R |] [(2.5.6.13)[link]]. In two dimensions with nonuniformly distributed projections, when the analytical form of the distribution of projections is not known, an appropriate approximation to [c(| R |,\Psi )] 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[link]),[\eqalignno{c( {| R |,\Psi } ) = | R |\,{\rm d}R\,{\rm d}\Psi \to c(R_j ,\Psi _i ) &= R_j {1 \over {2\pi }}{{\Psi _{i + 1} - \Psi _{i - 1} } \over 2} &\cr & = R_j {\Delta \Psi _i \over 4\pi }.&(2.5.6.29)} ]For a given set of angles the weights [c(R_j ,\Psi _i )] are easily computed (Fig. 2.5.6.5[link]).

[Figure 2.5.6.5]

Figure 2.5.6.5 | top | pdf |

Nonuniform distribution of projections. The projection weights for the reconstruction algorithms are chosen such that the backprojection integral becomes approximated by a Riemann sum and are equal to the angular length of an arc [\Delta \Psi _i ] (2.5.6.29)[link]. In Fourier space, projections of an object with real-space radius D form rectangles with width 1/D. In the exact filter backprojection reconstruction method, the weights are derived based on the amount the overlap between projections (2.5.6.29)[link].

In three dimensions, the weighting (2.5.6.29)[link] is applicable in a single-axis tilt data-collection 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)[link] 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)[link] 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[link].

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[link]), who suggested for a 2D case of non­uniformly distributed projections a heuristic weighting function,[c(R_j ,\Psi _i ) = {{R_j } \over {\textstyle\sum_l {\exp \left[ { - {\rm constant}( {| {\Psi _i - \Psi _l } |\bmod \pi } )^2 } \right]} }}.\eqno(2.5.6.30)]The weighting function (2.5.6.30)[link] 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)[link], which we consider optimal.

Radermacher et al. (1986[link]) 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)[link], a finite length of backprojection (Fig. 2.5.6.3[link]). Such a `truncated' backprojection is[\hat b_i ( {\bf{r}} ) = \varphi _{2_ i} ( {{\bf x}_\boldtau } ) * l( {\bf{r}} ),\quad \boldtau \perp {\bf x}\eqno(2.5.6.31)]with projection directions [\boldtau( {\theta ,\psi } )_i ] and[\eqalign{l( {\bf{r}} ) &= \delta ( {{\bf x}_\boldtau } )t( {z_\boldtau } ), \cr t( {z_\boldtau } ) &= \left\{ \matrix{1\hfill & -(D/2)\le z\le (D/2),\hfill\cr 0\hfill &{\rm otherwise},\hfill}\right.}\eqno(2.5.6.32)]where D is the diameter of the object or the length of the backprojection l. By taking the Fourier transform of (2.5.6.31)[link] and using the central section theorem (2.5.6.8)[link], we obtain a 3D Fourier transform of the backprojected projection,[\Phi _{3_i } ( {u,v,w} ) = \Phi _{2_i } ( {{\bf{u}}_\boldtau } )D\,{\rm sinc}( {D\pi w_\boldtau } ),\eqno(2.5.6.33)]and the 3D reconstruction is obtained by the inverse Fourier transform of the sum of contributions given by (2.5.6.33)[link],[{\scr F}[ {\hat b( {\bf{r}} )} ] = \textstyle\sum\limits_i \Phi _{2_i } ( {{\bf{u}}_\boldtau } )D\,{\rm sinc}( {D\pi w_\boldtau } ).\eqno(2.5.6.34)][z_\boldtau ] and [w_\boldtau ] are variables in real and Fourier spaces, respectively, both extending in the direction of the projection direction [\boldtau( {\theta ,\psi } )_i ].

In this analysis, the transfer function of the backprojection algorithm is obtained by setting [\Phi _{2_i } ( {{\bf{u}}_\boldtau } ) = 1], 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:[\varphi _3 ( {\bf{r}} ) = {\scr F}[ {\hat b( {\bf{r}} )} ]c( {u,v,w} ),\eqno(2.5.6.35)]with the weighting function given by[c( {u,v,w} ) = {1 \over {\textstyle\sum_i {D\,{\rm sinc}} ( {D\pi w_\boldtau } )}}.\eqno(2.5.6.36)]The general weighting function (2.5.6.36)[link] is consistent with analytical solutions (2.5.6.11)[link] and (2.5.6.14)[link], as it can be shown that by assuming infinite support [( {D \to \infty } )] and continuous and uniform distribution of projection directions, in two dimensions one obtains [c( {u,v} ) = ( {u^2 + v^2 } )^{-1/2}] (Radermacher, 2000[link]).

The derivation of (2.5.6.36)[link] is based on the analysis of continuous functions and its direct application to discrete data results in reconstruction artifacts; therefore, Radermacher (1992[link]) proposed attenuating the sinc functions in (2.5.6.36)[link] 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)[link]. The general weighted backprojection reconstruction is implemented in SPIDER (using exponent-based weighting functions) (Frank et al., 1996[link]), Xmipp (Sorzano et al., 2004[link]) and SPARX (Hohn et al., 2007[link]).

Harauz & van Heel (1986[link]) 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[link]), as follows from the central section theorem (2.5.6.8)[link]. 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 is[o_{il} (R) = T[{DR\sin ( {\Psi _i - \Psi _l } )} ],\eqno(2.5.6.37)]where 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)[link] are restricted such that [0 \le \Psi _i - \Psi _l \le (\pi/2)]. In this formulation, the overlap is limited to the maximum frequency,[R_{il}^{\max } = {2 \over {D\sin ( {\Psi _i - \Psi _l } )}}.\eqno(2.5.6.38)]Thus, the overlap function becomes[\let\normalbaselines\relax\openup 5pt o_{il} (R) =\left \{ \matrix{1-(R/R_{il}^{\rm max})\hfill & 0 \le R \,\lt\, u_{il}^{\rm max}\hfill\cr 0\hfill &R\,\gt\, u_{il}^{\rm max}\hfill}.\right.\eqno(2.5.6.39)]In effect, the weighting function, called by the authors an `exact filter', is[c(R,\Psi _i ) = {1 \over {1 + \textstyle\sum_{l, l \ne i } {o_{il} ( R )} }}.\eqno(2.5.6.40)]

The weighting function (2.5.6.40)[link] 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[link]). The method is conceptually simple and computationally efficient. However, (2.5.6.40)[link] does not approximate the correct weighting well for a uniform distribution of projections [i.e., it should yield [c(R,\Psi _i ) = R]]. This, as can be seen by integrating (2.5.6.39)[link] over the whole angular range, is not the case. The exact filter backprojection reconstruction is implemented in the IMAGIC (van Heel et al., 1996[link]), SPIDER (Frank et al., 1996[link]) and SPARX (Hohn et al., 2007[link]) packages.

The 3D reconstruction methods based on filtered backprojection are commonly used in single-particle 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)[link], (2.5.6.36)[link] and (2.5.6.38)[link]], so the performance of the reconstruction algorithm can be optimized for a particular data-collection geometry by changing the value of D (Paul et al., 2004[link]). 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)[link], (2.5.6.36)[link] and (2.5.6.40)[link] remain approximations of the correct weighting function (2.5.6.29)[link].

2.5.6.6. Direct Fourier inversion

| top | pdf |

Direct Fourier methods are based on the central section theorem (2.5.6.8)[link]. A set of the 2D Fourier transforms of projections yields an approximation to [\Phi _3] on a nonuniform 3D grid, and a subsequent numerical 3D inverse Fourier transform gives an approximation to [\varphi _3 ]. 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 [\Phi _3 ] on a nonuniform grid. In effect, the 3D inverse FFT is not applicable and an additional step of recovering samples of [\Phi _3 ] on a uniform grid from the available samples on a nonuniform grid is necessary.

One possibility is to resample the nonuniformly sampled version of [\Phi _3 ] 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[link]). 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[link]). 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[link]; Schomberg & Timmer, 1995[link]). The gridding-based direct Fourier reconstruction algorithm (GDFR) (Penczek et al., 2004[link]) was developed specifically for single-particle reconstruction. It comprises three steps:

  • (1) The first step, known as `gridding', involves calculating for the Fourier transform of each projection the convolution[\textstyle\sum\limits_i {{\scr F}[ w ] * ( {c{\scr F}[ {\varphi _{2_{i}} } ]} )}, \eqno(2.5.6.41)]where c are `gridding weights' designed to compensate for the nonuniform distribution of the grid points and [{\scr F}[ w ]] is an appropriately chosen convolution kernel. After processing all projections, this step yields samples of [{\scr F}[ w ] * {\scr F}[ {\varphi _3 } ]] on a Cartesian grid.

  • (2) 3D inverse FFT is used to compute[w\varphi _3 = {\scr F}^{ - 1}\left [ {{\scr F}[ w ] * {\scr F}[ {\varphi _3 } ]} \right]\eqno(2.5.6.42)]on a Cartesian grid.

  • (3) The weights are removed using division:[\varphi _3 =w\varphi_3/w.\eqno(2.5.6.43)]

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 [{\scr F}[ w ]] whose support in Fourier space is `small'. In order to assure good computational efficiency of the convolution step (2.5.6.41)[link], in GDFR this support was set to six Fourier voxels. In addition, in order to prevent division by zeros in (2.5.6.43)[link], the weighting function [w = {\scr F}^{ - 1}\left [ {{\scr F}[ w ]} \right]] must be positive within the support of the reconstructed object. A recommendable window function is the separable, bell-shaped Kaiser–Bessel window (O'Sullivan, 1985[link]; Jackson et al., 1991[link]; Schomberg & Timmer, 1995[link]) (Fig. 2.5.6.6[link]a):[\let\normalbaselines\relax\openup 5pt {\scr F}[ {w_{{\rm{KB}}} ( {\bf{u}} )} ] = \left\{ \matrix{\displaystyle\prod_{v=1}^d{\displaystyle I_0\{\pi D\alpha s_v[1-(u_v/s_v)^2]^{1/2}\}\over 2s_v}, \hfill &{\bf{u}} \in [ { - {\bf{s}},{\bf{s}}} ], \hfill\cr 0,\hfill &{\bf{u}}\,\, \notin\,\, [ { - {\bf{s}},{\bf{s}}} ] ,\hfill}\right.\eqno(2.5.6.44)]where [I_0 ] is the zero-order modified Bessel function and α = 1.25 is a parameter. The weighting function associated with [{\scr F}[ {w_{{\rm{KB}}} } ]] is[w_{{\rm{KB}}} ( {\bf x} ) = \displaystyle\prod\limits_{\nu = 1}^d{\sinh\left[\pi D\alpha s_v\left(1-\{x_v/[\alpha(D/2)]^2\}\right)^{1/2}\right]\over \pi D\alpha s_v\left(1-\{x_v/[\alpha(D/2)]^2\}\right)^{1/2}}.\eqno(2.5.6.45)]Finally, the gridding weights c are chosen such that in discrete implementation of (2.5.6.41)[link] and (2.5.6.42)[link] we obtain a Riemann sum approximating the respective integral (Schomberg, 2002[link]). The Voronoi diagram (Okabe et al., 2000[link]) 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.7[link]a). In cryo-EM, 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 [{\scr F}[ {\varphi _3 } ]] 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[link]) the 2D reverse gridding method is used to compute the Fourier transform of each projection on a 2D polar grid. In this way, [{\scr F}[ {\varphi _3 } ]] 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.7[link]b). This significantly reduces the memory requirements and improves the computational efficiency of the method, particularly when a fast [O( {n\log n} )] algorithm for computing the Voronoi diagram on the unit sphere is employed (Renka, 1997[link]).

[Figure 2.5.6.6]

Figure 2.5.6.6 | top | pdf |

Kaiser–Bessel window function used in the gridding reconstruction algorithm GDFR. (a) In Fourier space, the window function is effectively zero outside the support of six Fourier pixels. (b) In real space, the zeros of the window function are beyond the radius of the reconstructed object.

[Figure 2.5.6.7]

Figure 2.5.6.7 | top | pdf |

Partitions of the sampling Fourier space using Voronoi diagrams in direct inversion algorithms. (a) Central section of a 3D Fourier volume with sampling points originating from 2D Fourier transforms of projections. Although 2D projections are sampled on a uniform (Cartesian) grid, the arbitrary rotations of projections in 3D space yields a nonuniform distribution of points in three dimensions. In effect, the 3D reconstruction by direct inversion using 3D FFT is not possible. (b) Voronoi diagram on a sphere in the GDFR algorithm. Using the reverse gridding method, the 2D Fourier transforms of projections are resampled onto 1D central lines using a constant angular step. In 3D Fourier space, they are located on central sections and their angular directions are evenly distributed on grand circles. However, since central sections have nonuniform distributions, the distribution of angular directions (sampling points on the unitary sphere) is also nonuniform and effectively random.

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)[link] 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)[link]–(3)[link][link] of the gridding method:

  • (1) the input image is divided by the weighting function [\varphi_2/w];

  • (2) the image is padded with zeros and 2D FFT is used to compute [{\scr F}[\varphi_2/w]];

  • (3) gridding is used to compute [{\scr F}[w]*{\scr F}[\varphi_2/w]] on an arbitrary nonuniform grid.

Note the reverse gridding does not require explicit gridding weights, as in the third step they are constant (Penczek et al., 2004[link]).

In GDFR, the third step of the reversed gridding results in a set of 1D Fourier central lines [{\scr F}[ {\varphi _{1l} } ]] calculated using a constant angular step. Clearly, upon inverse Fourier transform they amount to a Radon transform [\varphi _1 ( {u,\psi } )] of the projection and, upon repeating the process for all available projections, they yield a Radon transform [\varphi _1 ( u,{\boldtheta } )] of a 3D object [\varphi _3 ] (albeit non­uniformly 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[link]).

The results of a comparison of selected reconstruction algorithms are shown in Fig. 2.5.6.8[link]. 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[link]). 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 gridding-based 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[link]). 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 intermediate-frequency ranges. The general weighted backprojection with exponent-based weighting function (WBP1) performs quite well, although reproduction in the low-frequency 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 high-quality 3D reconstructions from nonuniformly distributed projections. SIRT yields a reconstruction that in the low- and intermediate-frequency 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[link](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 Fourier-space-based method used to generate the test data.] In terms of computational efficiency, gridding-based algorithms outperform weighted backprojection algorithms by a small factor while SIRT is approximately ten times slower (depending on the number of iterations used).

[Figure 2.5.6.8]

Figure 2.5.6.8 | top | pdf |

(a) Plot of correlation coefficients (CCs) calculated between Fourier transforms of the reconstructed structure and the original phantom as a function of the magnitude of the spatial frequency using five reconstruction algorithms. GDFR: gridding direct Fourier reconstruction with Voronoi weights; GD3D: gridding reconstruction with simplified gridding weights; WBP1: general weighted backprojection with exponent-based weighting function; WBP2: exact filter weighted backprojection; SIRT: simultaneous iterative reconstruction algorithm. Noise-free projection data were computed in Fourier space using the reverse gridding method. (b) Rescaled version of the low-frequency range of (a). Note the different scales in (a) and (b). The horizontal axis is scaled in absolute frequency units with 0.5 equal to the Nyquist frequency.

2.5.6.7. 3D reconstruction of symmetric objects

| top | pdf |

Many objects imaged by EM have symmetries; the two types that are often met are helical (phage tails, F-actin, microtubules, myosin thick filaments, and bacterial pili and flagella) and point-group 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 period­icity and scattering density [\varphi ( {r,\psi ,z} )] is defined by the Fourier–Bessel transform:[\eqalignno{ \Phi (R, \Psi, Z) &= \sum\limits_{-\infty}^{+ \infty} \exp \left[ in \left(\Psi + {\pi \over 2}\right)\right] \int\limits_{0}^{\infty} \int\limits_{0}^{2 \pi} \int\limits_{0}^{l} \varphi (r, \psi, z) &\cr &\quad \times J_{n} (2 \pi rR) \exp [-i(n\psi + 2 \pi zZ)]r \ \hbox{d}r \ \hbox{d}\psi \ \hbox{d}z &\cr &= \sum\limits_{n} G_{n} (R, Z) \exp \left[ in \left(\Psi + {\pi \over 2}\right)\right]. &(2.5.6.46)\cr}]The inverse transform has the form [\rho (r, \psi, z) = {\textstyle\sum\limits_{n} \int} g_{n} (r, Z) \exp (in\psi) \exp (2 \pi izZ) \ \hbox{d}Z, \eqno(2.5.6.47)]so that [g_{n}] and [G_{n}] are the mutual Bessel transforms [G_{n} (R, Z) = {\textstyle\int\limits_{0}^{\infty}} g_{n} (r,Z) \,J_{n} (2 \pi rR) 2 \pi r \ \hbox{d}r\eqno(2.5.6.48)]and[g_{n} (r, Z) = {\textstyle\int\limits_{0}^{\infty}} G_{n} (R, Z)\, J_{n} (2 \pi rR) 2 \pi R \ \hbox{d}R.\eqno(2.5.6.49)]

Owing to helical symmetry, (2.5.6.48)[link] and (2.5.6.49)[link] contain only those of the Bessel functions that satisfy the selection rule (Cochran et al., 1952[link]) [l = mp + (nq/ N), \eqno(2.5.6.50)]where N, q and p are the helix symmetry parameters, [m = 0, \pm 1, \pm 2, \ldots]. Each layer l is practically determined by the single function [J_{n}] 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)[link]. 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[link]), 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[link]). Further complications exist when the filament does not have a precisely defined helical symmetry, such as F-actin, which has a random variable twist (Egelman et al., 1982[link]). To address these problems, Egelman developed a real-space refinement method for the reconstruction of helical filaments that is capable of determining the helical symmetry of an unknown structure (Egelman, 2000[link]). In this approach, the Fourier–Bessel reconstruction algorithm is replaced by a general reconstruction algorithm discussed in Section 2.5.6.6[link] with the real-space projections supplied as input. Thus, the symmetry is not enforced within the reconstruction algorithm; instead, it is determined and imposed subsequently by real-space averaging.

The presence of point-group symmetry in the structure means that any projection yields as many symmetry-related (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[link]; Lawton & Prasad, 1996[link]; Liang et al., 2002[link]). 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)[link]–(2.5.6.49)[link]] (Crowther, Amos et al., 1970[link]; Crowther & Amos, 1972[link]).

The general reconstruction methods (algebraic, filtered backprojection, direct Fourier inversion) easily accommodate symmetries in the data. In major single-particle reconstruction packages, reconstruction programs are implemented such that the point-group 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[link]) an additional real-space 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 fan-beam 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 three-dimensional 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). Three-dimensional 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 three-dimensional 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 single-particle methods. Ultramicroscopy, 85, 225–234.
Egelman, E. H., Francis, N. & DeRosier, D. J. (1982). F-actin 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). Three-Dimensional 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). Three-dimensional reconstruction of icosahedral particles – the uncommon line. J. Struct. Biol. 116, 48–55.
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.
Grigorieff, N. (1998). Three-dimensional 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 X-rays. J. Math. Anal. Appl. 62, 1–23.
Harauz, G. & van Heel, M. (1986). Exact filters for general geometry three-dimensional 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 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.
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 sinogram-based 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 three-dimensional ray-driven 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 high-resolution 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 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.
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 three-dimensional 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). Three-dimensional reconstruction of single particles embedded in ice. Ultramicroscopy, 40, 33–53.
Penczek, P. A. (2002). Three-dimensional spectral signal-to-noise ratio for a class of reconstruction algorithms. J. Struct. Biol. 138, 34–46.
Penczek, P. A., Renka, R. & Schomberg, H. (2004). Gridding-based direct Fourier inversion of the three-dimensional ray transform. J. Opt. Soc. Am. A, 21, 499–509.
Penczek, P. A., Zhu, J. & Frank, J. (1996). A common-lines 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). Three-dimensional reconstruction with contrast transfer function compensation from defocus series. Scan. Microsc. Suppl. 11, 1–10.
Radermacher, M. (1992). Weighted back-projection methods. In Electron Tomography, edited by J. Frank, pp. 91–115. New York: Plenum.
Radermacher, M. (1994). Three-dimensional reconstruction from random projections: orientational alignment via Radon transforms. Ultramicroscopy, 53, 121–136.
Radermacher, M. (2000). Three-dimensional 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 3-D 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 gridding-based 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., 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.
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). 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