Three-dimensional alignment of density maps in cryo-electron microscopy

Abstract

A common task in cryo-electron microscopy data processing is to compare three-dimensional density maps of macromolecules. In this paper, we propose an algorithm for aligning three-dimensional density maps, which exploits common lines between projection images of the maps. The algorithm is fully automatic and handles rotations, reflections (handedness), and translations between the maps. In addition, the algorithm is applicable to any type of molecular symmetry without requiring any information regarding the symmetry of the maps. We evaluate our alignment algorithm on publicly available density maps, demonstrating its accuracy and efficiency. The algorithm is available at https://github.com/ShkolniskyLab/emalign.

Impact Statement

This paper describes a fast algorithm for the rotational and translational alignment of three-dimensional density maps. Such alignment is an essential step in cryo-electron microscopy data processing. The algorithm is available at https://github.com/ShkolniskyLab/emalign.

1. Introduction

Single particle cryo-electron microscopy (cryo-EM) is a method to determine the three-dimensional structure of biological macromolecules from their two-dimensional projection images acquired by an electron microscope^{( 1 )}. In this method, a sample of identical copies of the investigated molecule is quickly frozen in a thin layer of ice, where each copy is frozen at an unknown random orientation. The frozen sample is imaged by an electron microscope, resulting in two-dimensional images, where each image is a tomographic projection of one of the randomly oriented copies in the ice layer. The goal of single particle cryo-EM is to determine the three-dimensional structure of the molecule from the acquired two-dimensional images. A common task in cryo-EM data processing is to compare two density maps of the same molecule. This is required, for example, to estimate the resolution of the maps, evaluating their Fourier shell correlation curve^{( 2 )}, or to analyze their different conformations. All these tasks require to first align two density maps, that is, to orient them in the same way in a common coordinate system. Due to the nature of the cryo-EM imaging process, the two density maps may differ not only in their three-dimensional orientation (i.e., their “rotation”) but may also have different handedness (namely, reflected relative to each other) and may be centered differently with respect to a common coordinate system.

In this paper, we propose an algorithm for aligning two density maps, which is fully automatic and can handle rotations, translations, and reflections between the maps. The algorithm requires as an input only the two density maps. In particular, it does not assume knowledge of any other information such as the symmetry of the maps.

Formally, let $ {\phi}_1:{\mathrm{\mathbb{R}}}^3\to \mathrm{\mathbb{R}} $ and $ {\phi}_2:{\mathrm{\mathbb{R}}}^3\to \mathrm{\mathbb{R}} $ be two volumes such that

(1) $$ {\phi}_2(r)={\phi}_1\left( Or-t\right), $$

where $ r={\left(x,y,z\right)}^T\in {\mathrm{\mathbb{R}}}^3 $ , $ O\in O(3) $ and $ t={\left(\Delta x,\Delta y,\Delta z\right)}^T\in {\mathrm{\mathbb{R}}}^3 $ ( $ O(3) $ is the group of all orthogonal transformations of the three-dimensional space, namely rotations and reflections). The alignment problem is to estimate $ O $ and $ t $ given $ {\phi}_1 $ and $ {\phi}_2 $ . The matrix $ O $ is known as the orientation parameter, and the vector $ t $ as the translation parameter. In practice, we only get samples of $ {\phi}_1 $ and $ {\phi}_2 $ , arranged as three-dimensional arrays of size $ n\times n\times n $ , where $ n $ is the resolution of sampling. In cryo-EM, $ {\phi}_1 $ and $ {\phi}_2 $ represent two reconstructions of the same underlying molecule that we would like to compare (such as two half maps from a refinement process). In principle, it is possible to approximate the solution to the alignment problem using exhaustive search, by generating a set of candidate pairs $ \left({O}_i,{t}_i\right) $ , where $ {O}_i\in O(3) $ and $ {t}_i\in {\mathrm{\mathbb{R}}}^3 $ , and finding the pair which “best aligns” $ {\phi}_1 $ to $ {\phi}_2 $ in some chosen metric. The purpose of the alignment algorithm presented in this paper is to estimate the optimal alignment parameters in a fast and accurate way.

The paper is organized as follows. In Section 2, we review existing alignment algorithms. In Section 3, we give a high level simplified description of our algorithm. A detailed description is then given in Section 4. This description relies on a method for aligning a single projection image against a volume, a procedure which is described in Section 5. In Section 6, we discuss implementation considerations of the algorithm and analyze its complexity. An optional procedure for refining the estimated alignment parameters is described in Section 7. In Section 8, we demonstrate numerically the properties and performance of our algorithm. Finally, in Section 9, we discuss the properties and advantages of our algorithm.

2. Existing Methods

There exist several methods for three-dimensional alignment of molecular volumes. The Chimera software^{( 3 )} offers a semi-automatic alignment method which requires the user to approximately align the volumes manually and then refines this alignment using an optimization procedure. This means that a sufficiently accurate initial approximation for the alignment is required. Achieving this initial approximate alignment manually naturally takes time and effort, yet it is crucial for the success of Chimera’s alignment algorithm. The alignment procedure implemented by Chimera maximizes the correlation or overlap function between the two volumes by using a steepest descent optimization. The iterations of this optimization stop after reaching convergence or after 2,000 steps.

Another alignment method is the projection-based volume alignment (PBVA) algorithm^{( 4 )}. This method aligns a target volume to a reference volume by aligning multiple projections of the reference volume to the target volume whose orientation is unknown. The PBVA algorithm is based on finding two identical projections, a projection $ {P}_1 $ from the reference volume and a projection $ {P}_2 $ from the target volume as follows. The reference volume is projected at some known Euler angles, resulting in a projection $ {P}_1 $ , and the matching projection $ {P}_2 $ is found by maximizing the cross-correlation function between $ {P}_1 $ and a set of projections representing the possible projections of the target volume. The cross-correlation function is of five parameters—three Euler angles and two translation parameters in the plane of the projection $ {P}_2 $ ^{( 5 )}. Finally, the rotation between the volumes is estimated from the relation between the Euler angles corresponding to the projections $ {P}_1 $ and $ {P}_2 $ . After estimating the rotation between the volumes, the translation between them is found using projection images from the target volume. The translation between the volumes is estimated by least-squares regression using the two translation parameters of each projection from the target volume, where a minimum of two projections is required for calculating the three-dimensional translation vector. Using multiple projection images to estimate the translation between the volumes makes the alignment more robust.

The Xmipp software package^{( 6 )} also offers a three-dimensional alignment algorithm. It is based on expanding the two volumes using spherical harmonics followed by computing the cross-correlation function between the two spherical harmonics expansions representing the volumes^{( 7 )}. The process of expanding a volume into spherical harmonics is called the spherical Fourier transform (SFT) of the volume, where like the fast Fourier transform (FFT) algorithm, there exists an efficient algorithm for calculating the SFT^{( 7 )}. The process of calculating the cross-correlation function between the two spherical harmonics expansions of the volumes and estimating the rotation between the two volumes is implemented by a fast rotational matching (FRM) algorithm^{( 8 )}. After estimating the rotation between the two volumes, the translation between them is found by using the phase correlation algorithm^{( 9 )}.

Finally, the EMAN2 software package^{( 10 )} offers 2 three-dimensional alignment algorithms. In the first algorithm (implemented by the program e2align3d, now mostly obsolete), the rotation between the volumes is estimated using an exhaustive search for the three Euler angles of the rotation. First, the algorithm generates a set of candidate Euler angles with large angular increments. Then, the algorithm iteratively decreases the angular increments in the set of candidates in order to refine the resolution of the angular search^{( 11 )}. A much faster tree-based algorithm is implemented in the program e2proc3d. This method performs three-dimensional rotational and translational alignment using a hierarchical method with gradually decreasing downsampling in Fourier space. In Section 8, we compare our algorithm to this latter algorithm, as well as to the FRM algorithm implemented by Xmipp.

3. Outline of the Approach

We are given two volumes $ {\phi}_1 $ and $ {\phi}_2 $ satisfying (1). For simplicity, we assume for now that the volumes have no symmetry and are related by rotation only (no translation nor reflection). We generate a projection image from $ {\phi}_2 $ , denoted $ P $ , corresponding to an orientation given by a rotation matrix $ R $ . Since $ {\phi}_1 $ and $ {\phi}_2 $ are the same volume up to rotation, we can orient $ P $ relative to $ {\phi}_1 $ , that is, we can find the rotation $ \tilde{R} $ such that projecting $ {\phi}_1 $ in the orientation determined by $ \tilde{R} $ results in the image $ P $ . As we show below, it holds that $ \tilde{R}= OR $ , where $ O $ is the transformation from (1). Since $ R $ and $ \tilde{R} $ are known, we can estimate $ O $ as $ O=\tilde{R}{R}^T $ .

In practice, it may be that $ \tilde{R} $ is not determined uniquely by $ P $ , as for example, a volume may have two very similar views even if it is not symmetric. Moreover, the volumes to align are discretized and sometimes noisy, which introduces inaccuracies into the estimation of $ O $ . Thus, to estimate $ O $ more robustly, instead of using a single image $ P $ , we generate from $ {\phi}_2 $ multiple images $ {P}_1,\dots, {P}_N $ with orientations $ {R}_1,\dots, {R}_N $ , align each $ {P}_i $ to $ {\phi}_1 $ as above, resulting in estimates for $ O $ given by $ {O}_i={\tilde{R}}_i{R}_i^T $ , and then estimate $ O $ from all $ {O}_i $ simultaneously by solving.

$$ O=\underset{R}{\mathrm{argmin}}\sum \limits_{i=1}^N{\left\Vert {O}_i-R\right\Vert}_F^2, $$

where $ {\left\Vert \cdot \right\Vert}_F $ is the Frobenius matrix norm. In Section 4, we give an explicit solution for the latter optimization problem.

The key of the above procedure is estimating the orientation of a projection image $ P $ of $ {\phi}_2 $ in the coordinate system of $ {\phi}_1 $ . This is done by inspecting a large enough set of candidate rotations and finding the rotation $ \tilde{R} $ for which the induced common lines between $ P $ (when assuming its orientation is $ \tilde{R} $ ) and a set of projections generated from $ {\phi}_1 $ best agree. As inspecting each candidate rotation involves only one-dimensional operations (even if the input volumes are centered differently), it is very fast and highly parallelizable. Thus, this somewhat brute force approach is applicable to very large sets of candidate rotations (several thousands, for accurate alignment) and still results in a fast algorithm. We discuss below the complexity and advantages of this approach. We summarize the outline of our approach in Figure 1 and describe it in detail in Sections 4 and 5.

Figure 1. Outline of the algorithm.

In the above approach, we assume that $ O $ is a rotation. However, $ {\phi}_1 $ and $ {\phi}_2 $ may have a different handedness, and so $ O $ may include a reflection. The above approach can obviously be used to resolve the handedness by aligning $ {\phi}_2 $ to $ {\phi}_1 $ and to a reflected copy of $ {\phi}_1 $ , and determining whether a reflection is needed using some quality score of the alignment parameters (e.g., the correlation between the aligned volumes). However, as we show below, in our method, there is no need to actually align $ {\phi}_2 $ to a reflected copy of $ {\phi}_1 $ , saving roughly half of the computations (those required to actually align $ {\phi}_2 $ to a reflected copy of $ {\phi}_1 $ ), as explained in Section 4.

We next explain in detail the various steps of our algorithm, including handling translations, reflections, and symmetry in the volumes.

4. Estimating the Alignment Parameters

Consider two volumes $ {\phi}_1 $ and $ {\phi}_2 $ , where one volume is a rotated copy of the other (assuming for now no reflection nor translation between the volumes), namely (see (1))

(2) $$ {\phi}_2(r)={\phi}_1(Or),\hskip1em r={\left(x,y,z\right)}^T\in {\mathrm{\mathbb{R}}}^3, $$

where $ O $ is an unknown rotation matrix. Our goal is to find an estimate for $ O $ .

In case where $ {\phi}_1 $ and $ {\phi}_2 $ exhibit symmetry, the solution for $ O $ is not unique. To be concrete, we denote by $ SO(3) $ the group of all $ 3\times 3 $ rotation matrices. A group $ G\hskip0.35em \subseteq \hskip0.35em SO(3) $ is a symmetry group of a volume $ \phi $ , if for all $ g\in G $ it holds that

(3) $$ \hskip1em \phi (r)=\phi (gr),\hskip1em r={\left(x,y,z\right)}^T\in {\mathrm{\mathbb{R}}}^3. $$

In other words, a symmetry group of a volume is a group of rotations under which the volume is invariant (see^{( 12 )} for more details). If we denote the symmetry group of $ {\phi}_1 $ by $ {G}_1\hskip0.35em \subseteq \hskip0.35em SO(3) $ and define $ {r}^{\prime }= Or $ , then, from (2) and (3), we get for any symmetry element $ g\in {G}_1 $

(4) $$ {\phi}_2(r)={\phi}_1(Or)={\phi}_1\left({r}^{\prime}\right)={\phi}_1\left(g{r}^{\prime}\right)={\phi}_1(gOr). $$

Comparing the latter with (2), we conclude that the solution for $ O $ is not unique, and we thus replace the goal of finding $ O $ by finding any $ gO $ for some arbitrary element $ g\in {G}_1 $ of the symmetry group.

Note that we assume that $ O $ is a rotation, namely that $ {\phi}_1 $ and $ {\phi}_2 $ are related by rotation without reflection. The case where $ O $ is a reflection will be considered below. Let $ P $ be a projection image generated from $ {\phi}_2 $ using a rotation $ R $ , that is,

(5) $$ P(x,y)={\int}_{-\mathrm{\infty}}^{\mathrm{\infty}}{\phi}_2(Rr)dz={\int}_{-\mathrm{\infty}}^{\mathrm{\infty}}{\phi}_2({xR}^{(1)}+{yR}^{(2)}+{zR}^{(3)})dz, $$

where $ {R}^{(1)} $ , $ {R}^{(2)} $ , $ {R}^{(3)} $ are the columns of the matrix $ R $ and $ r={\left(x,y,z\right)}^T $ . From (2), we have that

(6) $$ {\phi}_2(Rr)={\phi}_1(ORr). $$

Thus, using (5), we have

(7) $$ P\left(x,y\right)=\int_{-\infty}^{\infty }{\phi}_2(Rr) dz=\int_{-\infty}^{\infty }{\phi}_1(ORr) dz. $$

Equation (7) implies that if $ P $ has orientation $ R $ with respect to $ {\phi}_2 $ , then it has orientation $ OR $ with respect to $ {\phi}_1 $ . In Section 5, we describe how to estimate $ OR $ given $ P $ and $ {\phi}_1 $ , namely how to estimate a rotation $ \tilde{R} $ that satisfies $ \tilde{R}= OR $ . If the volume $ {\phi}_1 $ is symmetric with symmetry group $ {G}_1 $ , then (as shown above) the rotation $ OR $ is equivalent to the rotation $ gOR $ for any $ g\in {G}_1 $ , and moreover, the two rotations cannot be distinguished. Thus, we conclude that

$$ \tilde{R}= gOR $$

for some unknown $ g\in {G}_1 $ . Using the latter equation, we can estimate $ O $ as

(8) $$ O={g}^T\tilde{R}{R}^T. $$

Note that in the latter equation $ R $ is known, $ \tilde{R} $ can be estimated using the algorithm in Section 5 below, and $ g $ can be arbitrary. Thus, (8) provides a means for estimating $ O $ .

However, to estimate $ O $ more robustly, we use multiple projections generated from $ {\phi}_2 $ . Let $ {R}_1,\dots, {R}_N $ be random rotations, and let $ {P}_1,\dots, {P}_N $ be the corresponding projections generated from $ {\phi}_2 $ according to (5). Using the procedure described above, we estimate for each $ {P}_i $ a rotation $ {\tilde{R}}_i $ that satisfies $ {\tilde{R}}_i={g}_i{OR}_i $ for some unknown $ {g}_i\in {G}_1 $ . Thus, as in (8), we can estimate $ O $ using any $ i\in \left\{1,\dots, N\right\} $ by

(9) $$ O={g}_i^T{\tilde{R}}_i{R}_i^T. $$

Contrary to (8), if we want the right-hand side of (9) to result in the same $ O $ for all $ i=1,\dots, N $ , then $ {g}_i $ cannot be arbitrary. In order to estimate $ O $ , we therefore need to find $ {g}_i $ , $ i=1,\dots, N $ and combine all estimates for $ O $ given in (9) into a single estimate.

To that end, define

(10) $$ {X}_i={R}_i{\tilde{R}}_i^T,\hskip1em i=1,\dots, N, $$

and look at the matrix $ H $ of size $ 3N\times 3N $ whose $ \left(i,j\right) $ block of size $ 3\times 3 $ is given by (see (9) and (10))

(11) $$ {\displaystyle \begin{array}{c}{H}_{ij}={X_i}^T{X}_j={\tilde{R}}_i{R}_i^T{R}_j{\tilde{R}}_j^T\\ {}={\tilde{R}}_i{\left({g}_i^T{\tilde{R}}_i\right)}^T{OO}^T\left({g}_j^T{\tilde{R}}_j\right){\tilde{R}}_j^T={g}_i{g}_j^T.\end{array}} $$

By a direct calculation, we get that the matrix of size $ 3N\times 3 $

(12) $$ \tilde{g}=\left(\begin{array}{c}{g}_1W\\ {}\vdots \\ {}{g}_NW\end{array}\right), $$

where $ W $ is an arbitrary $ 3\times 3 $ orthogonal matrix (i.e., $ {WW}^T={W}^TW=I $ ), satisfies

(13) $$ H\tilde{g}=N\tilde{g}. $$

Equation (11) also shows that the matrix $ H $ is of rank 3, which together with (13) implies that $ \tilde{g} $ can be calculated by arranging the three leading eigenvectors $ {v}_1 $ , $ {v}_2 $ , $ {v}_3 $ of $ H $ in a matrix

(14) $$ V={\left(\begin{array}{ccc}\mid & \mid & \mid \\ {}{v}_1& {v}_2& {v}_3\\ {}\mid & \mid & \mid \end{array}\right)}_{3N\times 3}=\left(\begin{array}{c}{V}_1\\ {}\vdots \\ {}{V}_N\end{array}\right), $$

whose $ 3\times 3 $ blocks $ {V}_1,\dots, {V}_N $ are $ {g}_1W,\dots, {g}_NW $ , for some unknown arbitrary $ W $ (see^{( 13 )} for a detailed derivation). In practice, at this point, we replace each $ {g}_iW $ by its closest orthogonal transformation, as described in^{( 14 )}, to improve its accuracy in the presence of noise and discretization errors.

Next, in order to extract an estimate for $ {g}_1,\dots, {g}_N $ from (14) (i.e., to eliminate $ W $ from the estimates in $ \tilde{g} $ given by (12)), we multiply each $ {g}_iW $ by $ {\left({g}_1W\right)}^T $ , resulting in

(15) $$ {g}_{est}=\left(\begin{array}{c}{g}_{est_1}\\ {}\vdots \\ {}{g}_{est_N}\end{array}\right)=\left(\begin{array}{c}{g}_1{WW}^T{g}_1^T\\ {}\vdots \\ {}{g}_N{WW}^T{g}_1^T\end{array}\right)=\left(\begin{array}{c}{g}_1{g_1}^T\\ {}\vdots \\ {}{g}_N{g_1}^T\end{array}\right). $$

Thus, each $ {g}_{est_i} $ is a rotation, even if $ W $ is not. We define $ {O}_i={g}_{est_i}^T{\tilde{R}}_i{R}_i^T $ , and using (9), we get for $ i=1,\dots, N $

(16) $$ {O}_i={g}_{est_i}^T{\tilde{R}}_i{R}_i^T={g}_1{g}_i^T{\tilde{R}}_i{R}_i^T={g}_1O. $$

Thus, we have $ N $ estimates for $ {g}_1O $ . Equation (4) states that $ {\phi}_2(r)={\phi}_1(Or)={\phi}_1(gOr) $ for any symmetry element $ g\in {G}_1 $ . Therefore, estimating $ {g}_1O $ is equivalent to estimating $ O $ . In order to estimate $ {g}_1O $ simultaneously from all $ {O}_i $ , $ i=1,\dots, N $ , we search for the rotation $ {O}_{est}^{(1)} $ (the superscript will be explained shortly) that satisfies

(17) $$ {O}_{est}^{(1)}=\underset{R}{\mathrm{argmin}}\sum \limits_{i=1}^N{\left\Vert {O}_i-R\right\Vert}_F^2. $$

In other words, $ {O}_{est}^{(1)} $ is the “closest” to all the estimated rotations $ {O}_i $ in the least squares sense. To solve (17), let $ \tilde{O} $ be the $ 3\times 3 $ matrix

(18) $$ \tilde{O}=\frac{1}{N}\sum \limits_{i=1}^N{O}_i. $$

In^{( 15 )}, it is proven that the solution to the optimization problem in (17) is

(19) $$ {O}_{est}^{(1)}={\tilde{U}\tilde{V}}^T, $$

where $ \tilde{O}=\tilde{U}\tilde{\Sigma}{\tilde{V}}^T $ is the singular value decomposition (SVD) of $ \tilde{O} $ . The algorithm for estimating $ {O}_{est}^{(1)} $ given $ {\phi}_1 $ and $ {\phi}_2 $ , as described above, is presented in Algorithm 1.

Algorithm 1 Estimating $ {O}_{est}^{(1)}. $

Input: Volumes $ {\phi}_1,{\phi}_2. $

1: Generate random rotations $ {\left\{{R}_i\right\}}_{i=1}^N $

2: Generate from $ {\phi}_2 $ projections $ {\left\{{P}_i\right\}}_{i=1}^N $ corresponding to the rotations $ {\left\{{R}_i\right\}}_{i=1}^N $ ▷ Eq. (5)

3: Apply Algorithm 2 to each $ {P}_i $ and $ {\phi}_1 $ . Denote the resulting rotations by $ {\left\{{\tilde{R}}_i\right\}}_{i=1}^N $

4: for $ i=1 $ to $ N $ do

5: $ \hskip2em \mathrm{Calculate}\hskip0.35em {X}_i={R}_i{\tilde{R}}_i^T $ ▷Eq. (10)

6: Construct the matrix $ 3N\times 3N $

$$ H=\left(\begin{array}{cccc}{I}_{3\times 3}& {X}_1^T{X}_2& \cdots & {X}_1^T{X}_N\\ {}{X}_2^T{X}_1& {I}_{3\times 3}& \cdots & {X}_2^T{X}_N\\ {}\vdots & \vdots & \ddots & \\ {}{X}_N^T{X}_1& {X}_N^T{X}_2& \cdots & {I}_{3\times 3}\end{array}\right) $$

7: Find the three leading eigenvectors $ {v}_1,{v}_2,{v}_3 $ of $ H $

8: Set

$$ V={\left(\begin{array}{ccc}\mid & \mid & \mid \\ {}{v}_1& {v}_2& {v}_3\\ {}\mid & \mid & \mid \end{array}\right)}_{3N\times 3}=\left(\begin{array}{c}{V}_1\\ {}\vdots \\ {}{V}_N\end{array}\right) $$

▷ Eq. (14)

9: Compute

$$ {g}_{est}=\left(\begin{array}{c}{g}_{est_1}\\ {}\vdots \\ {}{g}_{est_N}\end{array}\right)=\left(\begin{array}{c}{V}_1{V}_1^T\\ {}\vdots \\ {}{V}_N{V}_1^T\\ {}\end{array}\right) $$

▷ Eq. (15)

10: for $ i=1 $ to $ N $ do

11: Calculate $ {O}_i={g}_{est_i}^T{\tilde{R}}_i{R}_i^T $ ▷ Eq. (16)

12: Calculate $ \tilde{O}=\frac{1}{N}{\sum}_{i=1}^N{O}_i $ ▷ Eq. (18)

13: Calculate $ {O}_{est}^{(1)}={\tilde{U}\tilde{V}}^T $ , where $ \tilde{O}=\tilde{U}\tilde{\Sigma}{\tilde{V}}^T $ is the SVD of $ \tilde{O} $ ▷ Eq. (19)

Output: $ {O}_{est}^{(1)} $ ▷ Estimated rotation

To handle the case where $ {\phi}_1 $ and $ {\phi}_2 $ have a different handedness (namely, related by reflection), we can of course apply Algorithm 1 to $ {\phi}_2 $ and a reflected copy of $ {\phi}_1 $ . However, this would roughly double the runtime of the estimation process, as the most time-consuming step in Algorithm 1 is step 3, whose complexity is $ O\left({n}^3\log n\right) $ operations for a volume of size $ n\times n\times n $ voxels (see Section 5).

Alternatively, it is possible to augment the above algorithm to handle reflections without doubling its runtime. In the case where there is a reflection between $ {\phi}_1 $ and $ {\phi}_2 $ , we need to replace the relation in (2) by the relation

(20) $$ {\phi}_2(r)={\phi}_1(OJr),\hskip1em J=\left(\begin{array}{ccc}1& 0& 0\\ {}0& 1& 0\\ {}0& 0& -1\end{array}\right). $$

Note that $ J $ in (20) is a reflection and that $ O $ is a rotation. Repeating the above derivation starting from (20) shows that to estimate $ O $ in this case, we can use the same $ {R}_i $ used above and the same estimates $ {\tilde{R}}_i $ obtained above (steps 1 and 3 of Algorithm 1), but this time we get that $ {O}_i={g}_i^T{\tilde{R}}_i{JR}_i^TJ $ (compare with (9)). Then, we set $ {X}_i={JR}_iJ{\tilde{R}}_i^T $ (compare with (10)) and proceed as above, resulting in an estimate $ {O}_{est}^{(2)} $ (compare with (17)), which corresponds to the optimal alignment parameters if $ {\phi}_1 $ and $ {\phi}_2 $ have opposite handedness. Once we have the two estimates $ {O}_{est}^{(1)} $ and $ {O}_{est}^{(2)} $ for the alignment parameters between $ {\phi}_1 $ and $ {\phi}_2 $ (without and with reflection), we estimate the translation corresponding to each of $ {O}_{est}^{(1)} $ and $ {O}_{est}^{(2)} $ using phase correlation^{( 9 )} (see Appendix C for details). This results in two sets of alignment parameters (rotation+translation). We then apply both sets of parameters to $ {\phi}_2 $ to align it with $ {\phi}_1 $ and pick the parameters for which $ {\phi}_2 $ after alignment has higher correlation with $ {\phi}_1 $ . We denote the estimated parameters by $ \left({O}_{est},{t}_{est}\right) $ .

5. Projection Alignment

It remains to show how to implement step 3 of Algorithm 1, that is, how to find the orientation of a projection $ P $ of $ {\phi}_2 $ with respect to the coordinate system of $ {\phi}_1 $ . Mathematically, we would like to solve the equation

(21) $$ P\left(x,y\right)=\int_{-\infty}^{\infty }{\phi}_1(Rr) dz,\hskip1em r={\left(x,y,z\right)}^T\in {\mathrm{\mathbb{R}}}^3 $$

for the unknown rotation $ R $ . A brute force approach of testing many candidate rotations in search for the $ R $ that (best) satisfies (21) is prohibitively expensive, as it requires to compute a projection of $ {\phi}_1 $ for each candidate rotation (this is essentially projection matching). We therefore take a different approach, whose cost for inspecting each candidate rotation is much lower (in fact requires $ O(n) $ operations to test each candidate rotation for a volume $ {\phi}_1 $ discretized into an array of size $ n\times n\times n $ ).

The idea is to generate several projection images from $ {\phi}_1 $ , and then, for each candidate rotation, to check the agreement of the common lines between $ P $ and the projections of $ {\phi}_1 $ , assuming the orientation of $ P $ is given by the candidate rotation. We estimate the rotation corresponding to $ P $ as the candidate rotation that results in the best agreement. We next formalize this method and then analyze its complexity.

We start by considering the case where there is no translation between $ P $ and $ {\phi}_1 $ , namely $ P $ and $ {\phi}_1 $ satisfy (21), and our goal is to estimate $ R $ given $ P $ and $ {\phi}_1 $ . We generate $ N $ projection images from $ {\phi}_1 $ ( $ N $ is typically small, see Section 8), denoted $ {P}_1^{(a)},\dots, {P}_N^{(a)} $ , with rotations $ {R}_1^{(a)},\dots, {R}_N^{(a)} $ chosen uniformly at random (note that we deliberately reuse the notation $ N $ used in Section 4, as explained below). We generate a set of candidate rotations $ S $ , over which we will search for the solution $ R $ of (21). The set $ S $ consists of a large number of approximately equally spaced rotations. See Appendix B for a detailed description of the construction of $ S $ .

We will assume for each candidate rotation $ Q\in S $ that $ P $ was generated using the rotation $ Q $ (i.e., we assume that $ R $ in (21) is equal to $ Q $ ), compute the mean correlation of the common lines between $ P $ and $ {P}_1^{(a)},\dots, {P}_N^{(a)} $ , and choose as an estimate for $ R $ the rotation $ Q $ for which the mean correlation is highest. Specifically, for each $ Q\in S $ and $ {R}_i^{(a)} $ , $ i=1,\dots, N $ , we compute the direction of the common line between $ P $ and $ {P}_i^{(a)} $ , given by the angles $ {\alpha}_i $ in $ P $ and $ {\beta}_i $ in $ {P}_i^{(a)} $ , as explained in Appendix A. The common line property^{( 16 )} states that if $ Q=R $ then

$$ \hat{P}\left(\xi \cos {\alpha}_i,\xi \sin {\alpha}_i\right)={\hat{P}}_i^{(a)}\left(\xi \cos {\beta}_i,\xi \sin {\beta}_i\right),\hskip1em \xi \in \mathrm{\mathbb{R}}, $$

where $ \hat{P} $ and $ {\hat{P}}_i^{(a)} $ are the Fourier transforms of $ P $ and $ {P}_i^{(a)} $ , respectively (see Appendix A for a review of common lines and their properties). We thus define

$$ {\displaystyle \begin{array}{l}{f}_i\left(Q,\xi \right)=\hat{P}\left(\xi \cos {\alpha}_i,\xi \sin {\alpha}_i\right),\\ {}{g}_i\left(Q,\xi \right)={\hat{P}}_i^{(a)}\left(\xi \cos {\beta}_i,\xi \sin {\beta}_i\right),\end{array}} $$

and the cost function

(22) $$ \rho (Q)=\frac{1}{N}\mathrm{\Re}\sum \limits_{i=1}^N\frac{\int_0^{\mathrm{\infty}}{\overline{f}}_i(Q,\xi ){g}_i(Q,\xi )\hskip0.1em d\xi }{{\Vert {f}_i\Vert}_{L^2}{\Vert {g}_i\Vert}_{L^2}}, $$

where $ {\overline{f}}_i $ denotes the complex conjugate of $ {f}_i $ . In other words, $ \rho (Q) $ measures how well the common lines induced by $ Q $ between $ P $ and $ {P}_1^{(a)},\dots, {P}_N^{(a)} $ agree. We then set our estimate for $ R $ to be

$$ {R}_{est}=\underset{Q\in S}{\mathrm{argmax}}\rho (Q). $$

We explore the appropriate value for $ N $ in Section 8.

We now extend the above scheme to the case where $ P $ is not centered with respect to $ {\phi}_1 $ , namely $ P $ is given by

(23) $$ P\left(x-\Delta x,y-\Delta y\right)=\int_{-\infty}^{\infty }{\phi}_1(Rr) dz,\hskip1em r={\left(x,y,z\right)}^T\in {\mathrm{\mathbb{R}}}^3, $$

for an unknown rotation $ R $ and an unknown translation $ \left(\Delta x,\Delta y\right) $ . The idea for estimating $ R $ is the same as before, except that the calculation of the common lines should take into account the unknown translation, as we describe next.

We denote the unshifted version of $ P $ by $ \tilde{P} $ , which is given by

(24) $$ \tilde{P}\left(x,y\right)=\int_{-\infty}^{\infty }{\phi}_1(Rr) dz,\hskip1em r={\left(x,y,z\right)}^T\in {\mathrm{\mathbb{R}}}^3 $$

(this is exactly (21), but we repeat it to clearly set up the notation). Then,

$$ P\left(x,y\right)=\tilde{P}\left(x+\Delta x,y+\Delta y\right). $$

Taking the Fourier transform of both sides of the latter equation, we get that^{( 1 )}

(25) $$ \hat{P}\left({\omega}_x,{\omega}_y\right)=\hat{\tilde{P}}\left({\omega}_x,{\omega}_y\right){e}^{i\left({\omega}_x\Delta x+{\omega}_y\Delta y\right)}. $$

Suppose that the common line between $ \tilde{P} $ and $ {P}_i^{(a)} $ is given by the angles $ {\alpha}_i $ in $ \tilde{P} $ and $ {\beta}_i $ in $ {P}_i^{(a)} $ (see Appendix A). By definition of the common line, it holds that

$$ \hat{\tilde{P}}\left(\xi \cos {\alpha}_i,\xi \sin {\alpha}_i\right)={\hat{P}}_i^{(a)}\left(\xi \cos {\beta}_i,\xi \sin {\beta}_i\right). $$

Using (25), we get that

$$ \hat{P}\left(\xi \cos {\alpha}_i,\xi \sin {\alpha}_i\right){e}^{- i\xi \Delta \xi }={\hat{P}}_i^{(a)}\left(\xi \cos {\beta}_i,\xi \sin {\beta}_i\right), $$

where $ \Delta \xi =\Delta x\cos {\alpha}_i+\Delta y\sin {\alpha}_i $ is the one-dimensional shift between the projections along their common line. We assume that this one-dimensional shift is bounded by some number $ d $ .

Thus, we need to modify our cost function (22) to take into account also the unknown (one-dimensional) phase $ {e}^{- i\xi \Delta \xi } $ . We therefore define (with a slight abuse of notation in reusing the previous notation for the cost function)

$$ {\displaystyle \begin{array}{c}{f}_i\left(Q,\Delta \xi, \xi \right)=\hat{P}\left(\xi \cos {\alpha}_i,\xi \sin {\alpha}_i\right){e}^{- i\xi \Delta \xi },\\ {}{g}_i\left(Q,\xi \right)={\hat{P}}_i^{(a)}\left(\xi \cos {\beta}_i,\xi \sin {\beta}_i\right)\end{array}} $$

and the cost function

(26) $$ \rho (Q,\Delta \xi )=\frac{1}{N}\mathrm{\Re}\sum \limits_{i=1}^N\frac{\int_0^{\mathrm{\infty}}{\overline{f}}_i(Q,\Delta \xi, \xi ){g}_i(Q,\xi )\hskip0.1em d\xi }{{\Vert {f}_i\Vert}_{L^2}{\Vert {g}_i\Vert}_{L^2}}, $$

and set our estimate for the solution $ R $ of (23) to be

(27) $$ {R}_{est}=\underset{\begin{array}{c}Q\in S,\Delta \xi \in \left[-d,d\right]\end{array}}{\mathrm{argmax}}\rho \left(Q,\Delta \xi \right). $$

The formula for the angles $ {\alpha}_i $ and $ {\beta}_i $ of the common line between $ P $ and $ {P}_i^{(a)} $ induced by the rotations $ Q\in S $ and $ {R}_i $ is given in Appendix A. Note that at this point we are only interested in $ {R}_{est} $ and not in the translation $ \left(\Delta x,\Delta y\right) $ in $ P $ , as the relative translation between $ {\phi}_1 $ and $ {\phi}_2 $ is efficiently determined using phase correlation (see^{( 9 )} and Appendix C) once we have determined their relative rotation. The algorithm for solving equation (23) is summarized in Algorithm 2.

Algorithm 2 Projection alignment.

Input: Projection $ P $ and volume $ {\phi}_1 $ satisfying (23).

1: Generate random rotations $ {R}_1,\dots, {R}_N $

2: Generate from $ {\phi}_1 $ projections $ {P}_1^{(a)},\dots, {P}_N^{(a)} $ corresponding to the rotations $ {R}_1,\dots, {R}_N $

3: Generate candidate rotations $ S $ ▷ Appendix B

4: Compute

$$ {R}_{est}=\underset{Q\in S,\Delta \xi \in \left[-d,d\right]}{\mathrm{argmax}}\rho \left(Q,\Delta \xi \right). $$

▷ Eqs. 26 and (27)

Output: $ {R}_{est} $ ▷ Estimated rotation

As mentioned above, we use the same $ N $ in Sections 4 and 5. While, in principle, the number of projections generated from $ {\phi}_2 $ in Section 4 can be different from the number of projections generated from $ {\phi}_1 $ in Section 5, due to the symmetric role of $ {\phi}_1 $ and $ {\phi}_2 $ in the alignment problem, there is no reason to consider different values.

6. Implementation and Complexity Analysis

Algorithms 1 and 2 are formulated in the continuous domain. Obviously, to implement them, we must explain how to apply them to volumes $ {\phi}_1 $ and $ {\phi}_2 $ given as three-dimensional arrays of size $ n\times n\times n $ . We now explain how to discretize each of the steps of Algorithms 4 and 5 and analyze their complexity. For simplicity, we use for the discrete quantities the same notation we have used for the continuous ones.

The only step in Algorithm 1 that needs to be discretized is step 2. This step is accurately discretized based on the Fourier projection slice theorem (A.1) using a non-equally spaced FFT^{( 17 , 18 )}, whose complexity is $ O\left({n}^3\log n\right) $ (for a fixed prescribed accuracy). The result of this step is a discrete projection image $ P $ given as a two-dimensional array of size $ n\times n $ pixels. The remaining steps of Algorithm 4 are already discrete, and since the value of $ N $ is small compared to $ n $ , their complexity is negligible.

We next analyze Algorithm 2. The input to this algorithm is a projection image $ P $ of size $ n\times n $ pixels, and a volume $ {\phi}_1 $ of size $ n\times n\times n $ voxels. The algorithm also uses the parameter $ N $ , but since it is a small constant, we ignore it in our complexity analysis. Step 1 of Algorithm 2 requires a constant number of operations. Step 2 is accurately implemented using a non-equally spaced FFT^{( 17 , 18 )}, whose complexity is $ O\left({n}^3\log n\right) $ (for a fixed prescribed accuracy). Step 3 is independent of the input volume, and moreover, the set $ S $ can be precomputed and stored. To implement step 4, we first discretize the interval of one-dimensional shifts $ \left[-d,d\right] $ in fixed steps of $ \Delta d $ pixels (say, 1 pixel). Specifically, we use the following shift candidates for the optimization in step 4:

$$ \Delta \xi \in \left\{-d+k\Delta d\;|\;k=0,\dots, \left\lfloor 2d/\Delta d\right\rfloor \right\}. $$

Then, for each $ Q\in S $ , we compute the angles $ {\alpha}_i $ and $ {\beta}_i $ (see Appendix A) and evaluate (26) for the pair $ \left(Q,\Delta \xi \right) $ by replacing the integral with a sum. If we store the polar Fourier transforms of all involved projection images $ P $ and $ {P}_1^{(a)},\dots, {P}_N^{(a)} $ (computed using the non-equally spaced FFT^{( 17 , 18 )}), each such evaluation amounts to accessing the rays in the polar Fourier transform corresponding to the angles $ {\alpha}_i $ and $ {\beta}_i $ , namely $ O(n) $ operations. Thus, the total number of operations required to implement step 4 of Algorithm 2 is $ \mid S\mid \times \left(\left\lfloor 2d/\Delta d\right\rfloor +1\right)\times n $ ( $ \mid S\mid $ is the number of elements in the set $ S $ ). Of course, all $ \mid S\mid \times \left(\left\lfloor 2d/\Delta d\right\rfloor +1\right) $ evaluations are independent and can be computed in parallel. Thus, the total complexity of Algorithm 2 is $ O\left({n}^3\log n\right) $ operations for step 2 and $ O(n) $ operations for testing each pair $ \left(Q,\Delta \xi \right) $ in step 4. Therefore, since the optimization in step 4 is very fast, it is practical to test even a very large set of candidate rotations $ S $ .

Finally, we note that in practice, to further speed up the algorithm, we first downsample the input volumes to size $ {n}_{ds} $ , align the two downsampled volumes, and apply the estimated alignment parameters to the original volumes. We demonstrate in Section 8 that this approach still results in a highly accurate alignment.

To understand the theoretical advantage of the above approach, we compare it to a brute force approach. In the brute force approach, we 1) scan over a large set of rotations and three-dimensional translations, 2) for each pair of a rotation and a translation, we transform one of the volumes according to this pair of parameters, and 3) choose the pair for which the correlation between the volumes after the transformation is maximal. Testing each pair of candidate parameters requires $ O\left({n}^3\right) $ operations (for rotating and translating one of the volumes, and for computing correlation), which amounts to a total of $ O\left({n}^3\times |S|\times {\left(2d/\Delta d\right)}^3\right) $ operations. In other words, testing each candidate rotation and translation is way more expensive than in our proposed method. In our approach, the expensive operation of complexity $ O\left({n}^3\log n\right) $ needs to be executed only once per each pair of inputs $ P $ and $ {\phi}_1 $ . Moreover, in our approach, the search over shifts is one-dimensional as opposed to the three-dimensional search required in the brute force approach.

7. Parameters’ Refinement

In this section, we describe an optional refinement procedure for improving the accuracy of the estimated parameters $ {O}_{est} $ and $ {t}_{est} $ obtained using the algorithm of Section 4.

We define the vector $ \Theta =\left(\psi, \vartheta, \varphi, {\Delta}_x,{\Delta}_y,{\Delta}_z\right) $ consisting of the six parameters required to describe the transformation between two volumes—three Euler angles ( $ \psi, \vartheta, \varphi $ ) describing their relative rotation and three parameters $ \left(\Delta x,\Delta y,\Delta z\right) $ describing their relative translation. We define the operator $ {T}_{\Theta}\left(\phi \right) $ , which applies the transformation parameters $ \Theta $ to the volume $ \phi $ (i.e., $ {T}_{\Theta} $ first rotates the volume and then translates it, according to the parameters in $ \Theta $ ). Next, for given volumes $ {\phi}_1 $ and $ {\phi}_2 $ , we denote their correlation by $ \rho \left({\phi}_1,{\phi}_2\right) $ . We are reusing the notation $ \rho $ from Section 5, since all occurrences of $ \rho $ in this paper correspond to a correlation coefficient whose evaluation formula is clear from its arguments. Finally, we define the objective function

(28) $$ c\left(\theta \right)=1-\rho \left({T}_{\Theta}\left({\phi}_1\right),{\phi}_2\right), $$

which vanishes for the parameters $ \Theta $ that align $ {\phi}_1 $ with $ {\phi}_2 $ .

To refine $ {O}_{est} $ and $ {t}_{est} $ of Section 4, we simply apply the BFGS algorithm^{( 19 )} to the objective function (28), with an initialization given by $ {O}_{est} $ and $ {t}_{est} $ .

8. Results

The alignment algorithm (with and without the optional refinement described in Section 7) was implemented in Python and is available online ¹, including the code that generates the figures of this section. A Matlab version of the algorithm is available as part of the ASPIRE software package^{( 20 )}.

As the algorithm uses two parameters—the downsampling $ {n}_{ds} $ (see Section 6) and the number of reference projections $ N $ (see Section 4)—we first examine how to appropriately set their values. Then, we examine the advantage of the refinement procedure proposed in Section 7. To show the benefits of our algorithm in practice, we then compare its performance to that of two other alignment algorithms—the alignment algorithm from the EMAN2 software package^{( 10 )} (implemented in the program e2proc3d) and the FRM algorithm implemented in the Xmipp software package^{( 6 )}. Finally, we examine the performance of the three algorithms using noisy input volumes.

We tested our algorithm on volumes from the electron microscopy data bank (EMDB)^{( 21 )} with different types of symmetries, whose properties are described in Table 1. All tests were executed on a dual Intel Xeon E5-2683 CPU (32 cores in total), with 768 GB of RAM running Linux. The memory required by the algorithm is of the order of the size of the input volumes. We used $ \mathrm{15,236} $ candidate rotations in Algorithm 2 (the size of the set $ S $ ), generated as described in Appendix B. This set of candidates is roughly equally spaced in the set of rotations $ SO(3) $ . While it is difficult to characterize the resolution of this set in terms of the resolution of each of the Euler angles, a rough calculation suggests that the resolution in each of the Euler angles is smaller than 5 degrees. We do not use rotations generated by a regular grid of Euler angles, as such a grid is less efficient than our grid, due to the nonuniform rotations generated by a regular grid of Euler angles. For example, discretizing each of the Euler angles to 5 degrees would result in 186,624 rotations, more than an order of magnitude larger than the number of rotations we use.

Table 1. Test volumes.

Note. Each volume is a three-dimensional array of size $ n\times n\times n $ , with $ n $ specified on the third column. The symmetry of each volume is given by the second column.

For each test, we generate a pair of volumes $ {\phi}_1 $ and $ {\phi}_2 $ related by a rotation matrix $ O $ and a translation vector $ t\in {\mathrm{\mathbb{R}}}^3 $ . The translation is chosen at random with magnitude up to 10% of the size of the volume. We denote the alignment parameters estimated by our algorithm by $ {O}_{est} $ and $ {t}_{est} $ . We evaluate the accuracy of our algorithm by calculating the difference between the rotations $ O $ and $ {O}_{est} $ . To that end, we first note that following (4), $ {O}_{est} $ is an estimate of $ gO $ for some arbitrary $ g\in {G}_1 $ , where $ {G}_1\hskip0.35em \subseteq \hskip0.35em SO(3) $ is the symmetry group of $ {\phi}_1 $ . In order to calculate the difference between $ O $ and $ {O}_{est} $ , we have to find the symmetry element $ g $ . In our tests, the symmetry group $ {G}_1 $ is known (see Table 1), and so we find $ g $ by solving

(29) $$ \underset{g\in {G}_1}{\mathrm{argmin}}{\left\Vert {O}_{est}- gO\right\Vert}_F, $$

followed by defining $ {O}_{est}^{\prime }={g}^T{O}_{est} $ . Next, the error in the estimated rotation $ {O}_{est}^{\prime } $ is calculated using the axis-angle representation of rotations as follows. The axis of the rotation $ O $ is defined to be the unit vector $ v\in {\mathrm{\mathbb{R}}}^3 $ that satisfies $ Ov=v $ , that is, $ v $ is an eigenvector of $ O $ corresponding to eigenvalue $ 1 $ . Similarly, we define the unit vector $ {v}^{\prime}\in {\mathrm{\mathbb{R}}}^3 $ to be the axis of the rotation $ {O}_{est}^{\prime } $ . Then, we calculate the angle between the axes of the rotations as

(30) $$ {e}_1={\cos}^{-1}\left({v}^T{v}^{\prime}\right). $$

The angle of rotation of the matrix $ O $ around its axis $ v $ is given by $ {\theta}_1={\cos}^{-1}\left(u\cdot Ou\right) $ , where $ u\in {\mathrm{\mathbb{R}}}^3 $ is a unit vector perpendicular to $ v $ . Similarly, we define $ {\theta}_2 $ to be the angle of rotation of the matrix $ {O}_{est}^{\prime } $ around its axis $ {v}^{\prime } $ . The error in the rotation angle is then defined as

(31) $$ {e}_2=\mid {\theta}_1-{\theta}_2\mid . $$

We start by investigating the appropriate value for the downsampling parameter $ {n}_{ds} $ (see Section 6). To that end, for each of the volumes in Table 1, we create its rotated and shifted copy and apply our algorithm with the downsampling parameter equal to 16, 32, 64, and 128 (namely, we actually align downsampled copies of the volumes and then apply the estimated parameters to the original volumes). The results are shown in Figure 2. For each value of downsampling, we show a bar plot that summarizes the results for all test volumes. Note that these results are without the refinement procedure of Section 7. To provide a more detailed information on the chosen downsampling value, we show in Figure 3 only the results for downsampling to sizes 64 and 128. Based on these results, we use a downsampling value of 64 in all subsequent tests. In particular, this value of downsampling results in an accurate initialization of the refinement procedure of Section 7, as shown in Figure 4. As of timing, we show in Figure 5 the timing, without and with refinement, for downsampling to sizes 64 and 128.

Figure 2. Downsampling parameter versus accuracy of the algorithm. The left figure corresponds to the error $ {e}_1 $ in the rotation axis (see (30)). The right figure corresponds to the error $ {e}_2 $ in the rotation angle (see (31)).

Figure 3. Downsampling parameter versus accuracy of the algorithm, focused on 64 and 128. See Figure 2 for more details.

Figure 4. Error without (left figure) and with (right figure) refinement for downsampling to size $ 64\times 64\times 64 $ . The error reported in the figure is either $ {e}_1 $ (30) or $ {e}_2 $ (31), as shown in the $ x $ -axis.

Figure 5. Timing of the alignment algorithm with downsampling to sizes 64 and 128. NR stands for “without refinement”; R stands for “with refinement.”

Next, we wish to determine the number of reference projections $ N $ to use in Algorithms 1 and 2. We set the downsampling parameter to 64 and measure the estimation error for different numbers of reference projections. The results are summarized in Figure 6. We also show the timing for different numbers of reference projections, without and with refinement, in Figure 7. Based on these results, we choose the number of reference projections to be 30, as a good compromise between accuracy and speed.

Figure 6. Error versus the number of reference projections $ N $ . The left and right figures show the error without and with the refinement procedure of Section 7, respectively. The error reported in this figure is the sum $ {e}_1+{e}_2 $ given in (30) and (31).

Figure 7. Time versus the number of reference projections.

Next, we compare the performance of our algorithm with that of EMAN2’s and Xmipp’s alignment algorithms. The accuracy and timing results are summarized in Tables 2 and 3, respectively. Finally, we demonstrate the performance of the different algorithms for noisy input volumes. To that end, we use as a reference volume EMD 2660^{( 22 )} from EMDB (of size $ 360\times 360\times 360 $ voxels) and create its rotated and translated copy. We add to the reference volume and its rotate/translated copy additive Gaussian noise with SNR (signal-to-noise ratio) ranging from 1 to 1/256. A central slice from the noisy reference volume at different levels of SNR is shown in Figure 8. The accuracy results of all algorithms for the various SNRs are shown in Table 4. The timings of the different algorithms are shown in Table 5.

Table 2. Accuracy comparison with EMAN2 and Xmipp.

Note. The errors reported in this table are the sum $ {e}_1+{e}_2 $ given in (30) and (31). Errors are given in degrees. For EMalign, (NR) corresponds to “without refinement” and (R) to “with refinement.” The two bottom rows show the mean and standard deviation of the error (in degrees) over all experiments.

Table 3. Timing comparison with EMAN2 and Xmipp (in seconds).

Note. For EMalign, (NR) corresponds to “without refinement” and (R) to “with refinement”. The column “size” is the side length of the input volumes.

Figure 8. Central slice of the noisy reference volume at different SNRs.

Table 4. Accuracy comparison for noisy input volumes at different SNRs.

Note. See Table 2 for more details.

Table 5. Timing comparison for noisy input volumes at different SNRs.

Note. All timings are given in seconds.

9. Discussion and Conclusions

In this paper, we proposed a fully automatic method for aligning three-dimensional volumes with respect to rotation, translation, and reflection. While the parameters of the algorithm can be tuned whenever needed, we showed that the default parameters work very well for a wide range of volumes of various symmetries. We also developed an auxiliary algorithm, which finds the orientation of a volume giving rise to a given projection image (Section 5). This algorithm may serve as a fast and highly accurate substitute to projection matching.

The core difference between our approach and other existing approaches is that our approach is based on common lines between projection images generated from the volumes. The advantage of this approach is that inspecting each candidate rotation is very fast, as it is based on one-dimensional operations on the common lines ( $ O(n) $ operations for volumes of size $ n\times n\times n $ ). We also note that our cost function (26) for identifying the optimal alignment is different than in other algorithms. While the typical cost function used by alignment algorithms is the correlation between the volumes, our cost function is the average correlation of the common lines between projection images of the volumes. These two cost functions are not equivalent, and while in our experiments we have not identified a scenario where one cost function is superior over the other, having tools that are based on different principles may prove beneficial in the future.

From the comparison of our algorithm with the alignment algorithms in EMAN2 and Xmipp, we conclude that our algorithm can be used in one of two modes. If we are interested in fast alignment with good accuracy (average error of 1.9 degrees of the rotation axis, and average error of 1.86 degrees of the in-plane rotation angle, with standard deviations of 1.25 degrees and 1.3 degrees, respectively), we can use our algorithm without the refinement procedure of Section 7. This is appropriate, for example, for visualization, as such an initial alignment is sufficient as an input for high resolution optimization-based alignment algorithms, such as the one in Chimera^{( 3 )}. In such a case, our algorithm is more than 3 times faster than EMAN2’s algorithm (even though our algorithm is implemented entirely in Python), and almost 40 times faster than Xmipp’s algorithm. If we are interested in very low alignment errors, the refinement procedure of Section 7 brings the average errors down to 0.25 degrees for the rotation axis and 0.28 degrees for the in-plane rotation angle (with standard deviations of 0.66 degrees and 0.63 degrees, respectively). In such a case, our algorithm is 80% faster than EMAN2’s and 15 times faster than Xmipp’s.

As of noise robustness, we see from Table 4 that our algorithm performs well down to SNR = 1/128. The algorithms in EMAN2 and Xmipp give very accurate results at even lower SNRs, but at the cost of a much higher running time. As a future research direction, there are several ways to improve the robustness of our algorithm to noise. First, as our algorithm is based on generating projection images of the volumes, it is possible to apply image denoising methods to the projection images. Then, it is possible to incorporate denoising into the common lines matching step (step 4 in Algorithm 2), for example, by denoising the common lines as one-dimensional signals or by incorporating frequency-dependent weights into (26). This is expected to significantly improve the robustness of our algorithm to noise while increasing its running time only slightly.