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.

2. Existing Methods

There exist several methods for three-dimensional alignment of molecular volumes. The Chimera software⁽Reference Pettersen, Goddard and Huang ^{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⁽Reference Yu, Snapp, Ruiz and Radermacher ^{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 $ ⁽Reference Radermacher ^{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⁽Reference de la Rosa-Trevín, Otón and Marabini ^{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⁽Reference Chen, Pfeffer, Hrabe, Schuller and Förster ^{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⁽Reference Chen, Pfeffer, Hrabe, Schuller and Förster ^{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⁽Reference Kovacs and Wriggers ^{8 )}. After estimating the rotation between the two volumes, the translation between them is found by using the phase correlation algorithm⁽Reference Kuglin and Hines ^{9 )}.

Finally, the EMAN2 software package⁽Reference Tang, Peng and Baldwin ^{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⁽Reference van Heel ^{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 (Reference Cucuringu, Singer and Cowburn13) 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⁽Reference Cucuringu, Singer and Cowburn ^{13 )} for a detailed derivation). In practice, at this point, we replace each $ {g}_iW $ by its closest orthogonal transformation, as described in⁽Reference Arun, Huang and Blostein ^{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 (Reference Arun, Huang and Blostein14) (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⁽Reference Singer and Shkolnisky ^{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⁽Reference Kuglin and Hines ^{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⁽Reference Shkolnisky and Singer ^{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⁽Reference Singer and Shkolnisky ^{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⁽Reference Kuglin and Hines ^{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⁽Reference Barnett, Magland and af Klinteberg ^{17 ,} Reference Barnett ^{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⁽Reference Barnett, Magland and af Klinteberg ^{17 ,} Reference Barnett ^{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⁽Reference Barnett, Magland and af Klinteberg ^{17 ,} Reference Barnett ^{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⁽Reference Press, Teukolsky, Vetterling and Flannery ^{19 )} to the objective function (28), with an initialization given by $ {O}_{est} $ and $ {t}_{est} $ .