1. Introduction
We study the problem of computing a compressed representation of the covariance matrix of 2D cryoEM images for the purpose of performing principal component analysis (PCA). More precisely, we consider an image formation model where the measurement is defined by
(1) 
where is a radial function, denotes convolution, is an image, and the noise term.
We are motivated by single particle cryoelectron microscopy (cryoEM) imaging, which is an important technique for determining the 3D structure of macromolecules. In particular, the single particle reconstruction (SPR) problem asks to recover the 3D structure of a macromolecule from noisy 2D images of its tomographic projections along unknown viewing angles. In cryoEM, the mathematical model is a special case of (1) and is of the form
(2) 
where are 3D spatial coordinates, is the point spread function, is the electrostatic potential of a molecule, is a 3D rotation, and
the noise term. In computational microscopy, it is typical to work with the Fourier transform of the point spread function, which is known as the contrast transfer function (CTF). In the simplest case, each measurement could correspond to a single fixed molecule potential function
; however, in general, we may assume that eachcould be a random variable representing a mixture of molecules, conformational heterogeneity, cases where the images are not perfectly centered, or other measurement imperfections
[20, 22].In general, each measurement can be associated with a different point spread function; however, in practice, a group of measurements, called a defocus group, can share a common point spread function. We assume that the measurements are grouped into defocus groups. Given and for , our goal is to estimate the 2D covariance function of the images
(3) 
where is a random variable from the same distribution as the images , and . We assume that the distribution of the images is invariant to inplane rotations (which is typically the case in cryoEM).
In cryoEM, the random variable is of the form , where the random variable is an unknown viewing angle, and the random variable is a molecule potential. There is generally no physical reason for a molecule to prefer one inplane rotation to another so distributions of random variables of this form are generally invariant to inplane rotations. In the field of cryoEM processing, the covariance function is simply referred to as the 2D covariance.
1.1. Motivation
The 2Dcovariance is an essential component of a number of computational techniques in cryoEM; we survey a few of these below.
First, we are motivated by PCA, which is a ubiquitous technique in statistics, data science, and computational mathematics and has applications to dimensionality reduction, denoising, visualization, among others. The principal components (that is, the top eigenvectors of the digitized covariance matrix) have a number of uses in the computational cryoEM pipeline. The subspace corresponding to the top eigenvectors of the covariance matrix identifies salient features of the dataset which enables, for instance, improved methods for image classification and visualization, such as Multivariate Statistical Analysis
[27, 26, 28, 29]. These techniques improve computational speed, since clustering becomes computationally easier in a space of reduced dimension, as well as accuracy, since dimensionality reduction by PCA amplifies the effective signaltonoise ratio because many coordinates for which noise dominates the signal are eliminated
[32, Table 3].Second, the covariance matrix has applications in the method of moments, a classical statistical inference method, applied to cryoEM
[13]. In this method, the 2D covariance is used to compute the similarly defined autocorrelation function of the underlying 3D structure. Under further assumptions such as sufficient nonuniformity of the distribution of the viewing angles [24] or sufficient sparsity of the molecular structure [2], this autocorrelation function determines the 3D density map either up to a finite list of possible structures or uniquely, respectively. This has been further developed into principled methods for ab initio estimation of cryoEM structures [24, 2], with reduced risk of userinduced model bias in the initial model. Alternatively, when additional information is available, for instance one [11] or two [19] noiseless projection images, or the 3D structure of a related, homologous structure [3, 5], the 3D density map is uniquely determined by the autocorrelation, without requiring any structural assumptions.Third, the covariance matrix has applications to denoising and CTFcorrecting projection images. Covariance Wiener Filtering (CWF) [4] is an approach which uses the classic Wiener filtering framework with the estimated covariance matrix to solve the image deconvolution and denoising problem. The technique represents images in a lower dimensional subspace that is formed from PCA using the estimated covariance matrix. The method then applies Wiener filtering to correct the CTFs and denoise the images in this reduced subspace.
Compared to the standard PCA problem, the cryoEM setting exhibits further computational challenges, since the estimation method also has to account for convolution with the point spread function, which destroys information of the resulting convolved function; see §2.2 for a more precise statement. On the other hand, the problem has additional symmetries making fast algorithms possible. In this paper, we present a new fast algorithm for estimating the covariance matrix that improves upon past approaches (especially when there are a large number of defocus groups) in terms of time and space complexity.
1.2. Main contribution
The main contribution of this paper is a new computational method for estimating the covariance (3) from measurements of the form (1) encoded by digitized images. The presented fast method has time complexity independent of the number of defocus groups. This is in contrast to past methods, where this complexity scales poorly with and involves operations [4], where is the number of iterations needed in a Conjugate Gradient step. Many modern cryoEM experimental datasets fall into the computationally challenging regime where scales with .
Our fast method hinges on a new fast and accurate method for expanding images into the FourierBessel basis, which provides a convenient way to handle convolution of radial functions (such as point spread functions) with images: namely, convolution with radial functions can be expanded as a diagonal operator operating on the basis coefficients [21].
The FourierBessel basis functions are harmonics
on the disk: the standing waves associated with the resonant frequencies of a disk shaped drum with a fixed boundary. More precisely, the harmonics on the disk are eigenfunctions of the Laplacian on the unit disk that satisfy Dirichlet boundary conditions. In computational mathematics, this basis is referred to as the FourierBessel basis, since the basis functions can be expressed as a product of a Bessel function and a complex exponential; see (
5) for a definition.Because of this simple structure, the covariance matrix of clean images can be estimated by a simple closedform solution, without using the (computationally expensive) conjugate gradient method from previous approaches. Simultaneously, the covariance matrix retains its block diagonal structure, meaning that its diagonal blocks can be estimated separately and independently, which altogether makes PCA fast.
We present numerical results of covariance estimation on synthetic and experimental data. Additionally, we show how the estimated covariance matrix can be used to denoise images using CWF, and perform PCA to visualize eigenimages from experimental data. Code implementing the method is publicly available online^{1}^{1}1Code is available at https://github.com/yunpengshi/fastcryoEMPCA. Moreover, our approach has the potential to generalize to settings beyond cryoEM, where PCA is used for signals estimated under more general group actions [1].
1.3. Organization
2. Methodology
2.1. Notation
For two matrices and , we denote their Hadamard (or entrywise) product by , the Hadamard division of and by and the th Hadamard power of by . These operations are defined elementwise by
(4) 
respectively. If is an
dimensional vector, then
denotes the matrix with along its diagonal, i.e., , and zeros elsewhere. If is a radial function, we write to mean that can be expressed as a function only of the magnitude of .2.2. Technical details
We make the following assumptions

[label=(A0)]

We assume that the point spread functions in the model (1) are radial functions; this implies that their Fourier transform (the CTFs) are also radial. In systems where astigmatism is present and the point spread function deviates slightly from a radial function, our approach can be used as an initial approximation that could be refined using the Conjugate Gradient method.

We assume that the underlying images in the model (1) are i.i.d. random variables whose distribution is invariant to inplane rotations.
Informally speaking, assumption 3 can be interpreted as saying that for any pair of frequencies and there is a point spread function whose Fourier transform does not vanish at or . Assumption 1 implies that is a radial function. Figure 1 shows the values of for each pair of radial frequency in log scale, for 1081 distinct CTF images of size , whose defocus values range from 0.81 to 3.87 , for an experimental dataset; see §4 for more details.
This assumption is much weaker than assuming that for all that does not vanish at any frequency. In the latter case, we could just use to invert each equation to get access to the underlying functions . If we had direct access to the underlying functions , then we could approximate the covariance function by the sample covariance
where is the sample mean. Indeed,
by the law of large numbers. To clarify why it is useful for
not to vanish, note that in the Fourier domain our measurement model can be expressed aswhere , and denote the Fourier transforms of and , respectively. Since taking the Fourier transform changes convolution to pointwise multiplication, if each Fourier transform for some , then we could estimate by first estimating from for , and then using the sample covariance matrix. However, in practice, the CTF is approximately a radial function with many zerocrossings, which means that multiplication by destroys information in the corresponding frequencies, making the restoration from a single image illposed.
2.3. FourierBessel
The main ingredient in our fast covariance estimation is a fast transform into a convenient and computationally advantageous basis, known as the FourierBessel basis (which consists of the harmonics on the disk: eigenfunctions of the Laplacian on the disk that obey Dirichlet boundary conditions). This specific choice of basis has a number of beneficial properties:

[label=()]

it is orthonormal,

it is ordered by frequency,

it is steerable i.e., images can be rotated by applying a diagonal transform to the basis coefficients,

it is easy to convolve with radial functions i.e., images can be convolved with radial functions by applying a diagonal transform to the basis coefficients.
These properties have made the FourierBessel basis a natural choice in a number of imaging applications [25, 4, 23, 32, 31] and will be central to the development of our fast covariance estimation method. In polar coordinates in the unit disk the FourierBessel basis functions are defined by
(5) 
where is a normalization constant, is the th order Bessel function of the first kind [6, §10.2], and is its th smallest positive root; the indices run over .
Recent work [21] has devised a new fast algorithm to expand images into FourierBessel basis functions. Informally speaking, given basis coefficients, the algorithm can evaluate the function on an grid in operations; the adjoint can be computed in the same number of operations, which makes iterative methods fast.
2.4. Key property of FourierBessel basis
A key property of the FourierBessel basis is that convolution with radial functions are diagonal transformations in any truncated basis expansion. More precisely, the following result holds [21, Lemma 2.3]: suppose that for some index set , and is a radial function. Then,
(6) 
where denotes orthogonal projection onto the span of , is the Fourier transform of , and is the th positive root of .
We emphasize that the weights of the diagonal transform in (6) are not the coefficients of in the disk harmonic expansion. Indeed, since is radial, it follows from (5) that the coefficients of in the basis satisfy when . Computing the weights from the coefficients of in the basis, would involve computing weighted sums of the Fourier transforms of the basis functions: , for some index set .
As an alternative to the disk harmonic basis expansion, one can consider simply taking the Fouriertransform of (1), which also leads to a diagonal representation of the convolution operator. However, the discrete Fourier transform does not have the steerability property, which is essential for the covariance estimation. Another attempt could be to use the polar Fourier transform. However, this representation is not invariant to arbitrary inplane rotations, but only to finitely many rotations as determined by the discretization spacing of the grid. These expansions are therefore unsuitable for the goal of this article and we instead use expansions into the FourierBessel basis, although other steerable bases could be considered [14, 15]. Table 1 summarizes the considerations that make the FourierBessel basis a natural choice of basis.
Basis  Orth.  Cont. steerable  Radial convolution diag.  Fast expansion^{2}^{2}2The basis expansion from Cartesian grid representation can be completed within operations up to log factors. 

Real  ✓  ✗  ✗  ✓ 
2D discrete Fourier  ✓  ✗  ✓  ✓ 
Polar Fourier  ✗  ✗ (discrete)  ✓  ✓ 
PSWF [14, 15]  ✓  ✓  ✗  ✗ 
FourierBessel [21]^{3}^{3}3The new expansion algorithm [21] improves the previous computational method [32] in terms of accuracy guarantees, computational complexity and the fact that it derives weights such that radial convolution is a diagonal operation.  ✓  ✓  ✓  ✓ 
2.5. Block diagonal structure
The steerable property of the FourierBessel basis implies that the 2Dcovariance will be block diagonal in this basis. A full representation of an matrix requires elements, but this is reduced to nonzero entries by the block diagonal structure. This block diagonal structure follows from the form of the basis functions and that the distribution of inplane rotations is assumed to be uniform. Indeed, suppose that For simplicity assume that has mean zero; subtracting the mean will only change the radial components, since the other components corresponding to nonvanishing angular frequencies have zero mean by merely averaging over all possible inplane rotations. By assumption 2, the covariance function in polar coordinates satisfies
(7) 
for all in . The covariance function in (3) can be expanded in a double FourierBessel basis expansion as
(8) 
where is the covariance matrix in the FourierBessel basis. Combining (7) and (8) and integrating over gives
(9) 
where is a Dirac function that is equal to if and zero otherwise, and we note that the second to last equality uses the fact that . Since the coefficients in the expansion (8) are unique, it follows from (9) that when . Hence, the covariance matrix has a block diagonal structure whose blocks consist of the indices for a given value of . In the following section, we show how these properties enable a fast method to estimate the covariance matrix.
2.6. Covariance Estimation
In the FourierBessel basis, (1) is written as
(10) 
where , , and are coefficient vectors of in the FourierBessel basis, respectively and is the vector encoding the convolution operator of §2.4, i.e., with components . The vectors are dimensional column vectors, where is the number of basis coefficients. We use this simple structure to obtain a closed form expression for the sample covariance matrix of the . We estimate this matrix by minimizing the discrepancy between the sample covariance and the population covariance; more precisely, the estimated covariance matrix is computed by solving the least squaresproblem
(11) 
where is defined by
(12) 
whose solution is
(13) 
The least squaressolution of (11) can be determined by the following system of linear equations
(14) 
where
(15) 
It follows that
(16) 
As discussed in §2.4, the covariance matrix is block diagonal in the FourierBessel basis. More precisely, the only nonzero elements of the matrix are those with . Therefore, the matrices in (15) and (16) need only be calculated for this subset of indices. Since there is a total of of these indices, this reduces the computational complexity compared to computing with the full matrices.
Note that the covariance matrix estimated from (16) may not be positive semidefinite due to subtraction of the term
. Therefore, when running the method in practice it is beneficial to use an eigenvalue shrinkage method. The computational cost of eigenvalue shrinkage for a matrix with our block structure is
, see for example [7, 4]. For completeness, we include this computational cost of eigenvalue shrinkage in our overall computational complexity. Informally speaking, the idea of eigenvalue shrinkage is to replace the term in (16) by , and then shrink and truncate the eigenvalues in a systematic way before Hadamard division, see [7] for details. The steps of the algorithm are summarized in Algorithm 1.The complexity of Step 1 of the algorithm is . The complexity of Step 1 is . The complexity of Step 1 is , since the number of basis coefficients is . The complexity of Step 1 is , since there are at most nonzero elements with indices , The complexity of Step 1 is , where the additional term comes from the computational complexity of eigenvalue shrinkage. Thus, the total complexity of the algorithm is .
2.7. Application to image denoising
3. Synthetic data results
We compare the timings of our fast method to previous approaches [4], for synthetic images generated from the 3D volume of SARSCoV2 (Omicron) spike complexes [10] (EMD32743), from the online EM data bank [16]. The original volume has size in each dimension, with pixel size 0.832 Å. We downsample the original volume to size , with
and 512, respectively, and show the computational times. To generate the synthetic noisy images, we first generate 10000 clean projection images of the 3D volume from random and uniformly distributed viewing directions. We next divide the set of clean images into a number of defocus groups, where the defocus values range from 1
to 4 . For all CTFs, we set the voltage as 300 and the spherical aberration as 2 . After convolving the images with their CTFs, we add colored noise with power spectral density up to a constant scale, where is the radial frequency. For both the previous method and ours, the CTFs and the noisy images are whitened before estimating the covariance. A few sample clean and noisy images are shown in Figure 5. All experiments were carried out on a machine with 72 cores. We note that our implementation does not take advantage of all cores since some of the packages we use are not fully optimized for parallel computing; improving this is a technical direction for future work.Figure 2 shows the time required to estimate the covariance matrices as a function of the number of defocus groups. Note that the runtime for the previous approach for the largest image and defocus group sizes is infeasibly large and that our fast method exhibits a speedup of up to three orders of magnitude.
Figure 3 shows the top six principal components estimated by our method and by traditional PCA using raw images, where , SNR and defocus groups, compared to traditional PCA on images. We use the sample covariance matrix of phaseflipped images in real space for the traditional PCA. For all methods, we use to denote times the eigenvalues of the eigenimages. The eigenimages from the traditional PCA look much noisier than ours, and contain artifacts that are due to imperfect CTF correction (see e.g. the circular artifacts in top 3 eigenimages in Figure 3). The eigenimages from the traditional PCA also fail to preserve the symmetries (see e.g. 6th eigenimage in Figure 3) that are present in our eigenimages, since they do not utilize the steerable basis and rotationaugmented images.
Figure 4 shows the quality of covariance estimation and image denoising when and using defocus groups. The quality of the covariance estimation is measured by the relative error in each angular frequency, which is defined as
(18) 
where and are respectively the th diagonal blocks of the clean and estimated covariance matrix, corresponding to all indices of the form . The clean covariance matrix was approximated by the sample covariance matrix of clean projection images. The performance of image denoising is measured by the Fourier ring correlation (FRC) between the clean and denoised images. Namely, for the th pair of clean and denoised images and , we first compute their Fourier coefficient vectors at radial frequency by the nonuniform FFT [8, 9, 18], where and . We then compute their averaged correlation for each :
where denotes the inner product between two complex vectors. The FRC is a realvalued quantity due to a symmetry property that arises since the images are realvalued.
In addition to the speedups of Figure 2, the left panel of Figure 4 demonstrates a slight increase in the estimation quality of our proposed algorithm, compared to the previous approach. This is possibly caused by improved accuracy in the FourierBessel basis approximation as well as improved accuracy by using the closedform expression (16) compared to the approximate conjugate gradient step of previous approaches. Similarly, on the right panel, the Fourier ring correlation between the denoised and clean images shows a slight performance increase. Figure 5 shows sample denoised images for different values of the signaltonoise ratio (SNR) where and . As a comparison, we show images denoised using the approach of this paper and the CWF method [4].
4. Experimental data results
We conclude by using our method on two experimental datasets obtained from the Electron Microscopy Public Image Archive [12], namely EMPIAR10028 [30], and EMPIAR10081 [17]. EMPIAR10028 is a dataset of the Plasmodium falciparum 80S ribosome bound to the antiprotozoan drug emetine whose 3D reconstruction is available in the EM data bank as EMD2660 [30]. The dataset contains 105247 motion corrected and picked particle images, from 1081 defocus groups, of size with 1.34 Å pixel size. EMPIAR10081 is a dataset of the human HCN1 hyperpolarizationactivated cyclic nucleotidegated ion channel, whose 3D reconstruction can be found in the EM data bank as EMD8511 [17]. The dataset contains 55870 motion corrected and picked particle images, from 53384 defocus groups, of size with 1.3 Å pixel size.
EMPIAR10028  

Methods  
oldCWF,  1415  598  27550  201  29764 
fastCWF,  768  5  2220  95  3088 
EMPIAR10081  

Methods  
oldCWF,  47  3369  46318  45  49779 
fastCWF,  35  8  93  27  163 
oldCWF,  363  NA  NA  NA  NA 
fastCWF,  169  34  1007  163  1373 
Computational times are shown in Table 2, showing a speedup of more than two orders of magnitude for the datasets with the largest number of distinct CTFs. Note that regular CWF on EMPIAR10081 encounters memory issues and cannot be run to completion, whereas our fast method runs seamlessly. In order to obtain a comparison, we therefore additionally downsample these images to where the original CWF can successfully run. On EMPIAR10028, we used all images for covariance estimation, and denoised the 2014 images from th, th, th, …, th defocus groups. On EMPIAR10081, we used all images for covariance estimation, and denoised the 502 images from th, th, th, …, th defocus groups. For both old and new methods, the covariance matrices are further refined to correct image contrast variations [25]. Sample visualization results are shown in Figures 6 and 7.
5. Discussion
Covariance estimation and PCA of cryoEM images are key ingredients in many classic cryoEM methods including multivariate statistical analysis [27, 26, 28, 29] and Kam’s method for abinitio modeling [13]. We propose a fast method to estimate the covariance matrix of noisy cryoEM images, and then illustrate its application to simultaneously correct for the CTFs and denoise the images. The approach relies on recent improvements to algorithms for expanding images in the FourierBessel basis [21], and has time complexity which is independent of the number of defocus groups. Our new approach is both significantly faster and more memoryefficient compared to the previous CWF method [4] and we apply our method to large experimental datasets with many distinct CTFs with speedups by factors up to more than two orders of magnitude. Our approach could potentially be extended to higherdimensional data and to the setting where images are distorted by CTFs which are not exactly radial, using either analytical correction terms or iterative numerical steps.
References
 [1] A. S. Bandeira, B. BlumSmith, J. Kileel, A. Perry, J. Weed, and A. S. Wein. Estimation under group actions: recovering orbits from invariants. arXiv preprint arXiv:1712.10163, 2017.
 [2] T. Bendory, Y. Khoo, J. Kileel, O. Mickelin, and A. Singer. Autocorrelation analysis for cryoEM with sparsity constraints: Improved sample complexity and projectionbased algorithms. arXiv preprint arXiv:2209.10531, 2022.
 [3] T. Bhamre, T. Zhang, and A. Singer. Orthogonal matrix retrieval in cryoelectron microscopy. In 2015 IEEE 12th International Symposium on Biomedical Imaging (ISBI), pages 1048–1052. IEEE, 2015.
 [4] T. Bhamre, T. Zhang, and A. Singer. Denoising and covariance estimation of single particle cryoEM images. Journal of structural biology, 195(1):72–81, 2016.
 [5] T. Bhamre, T. Zhang, and A. Singer. Anisotropic twicing for single particle reconstruction using autocorrelation analysis. arXiv preprint arXiv:1704.07969, 2017.
 [6] NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.1.5 of 20220315. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
 [7] D. L. Donoho, M. Gavish, and I. M. Johnstone. Optimal shrinkage of eigenvalues in the spiked covariance model. Annals of statistics, 46(4):1742, 2018.
 [8] A. Dutt and V. Rokhlin. Fast fourier transforms for nonequispaced data. SIAM Journal on Scientific computing, 14(6):1368–1393, 1993.
 [9] L. Greengard and J.Y. Lee. Accelerating the nonuniform fast Fourier transform. SIAM Review, 46(3):443–454, 2004.
 [10] H. Guo, Y. Gao, T. Li, T. Li, Y. Lu, L. Zheng, Y. Liu, T. Yang, F. Luo, S. Song, et al. Structures of omicron spike complexes and implications for neutralizing antibody development. Cell reports, 39(5):110770, 2022.
 [11] S. Huang, M. Zehni, I. Dokmanić, and Z. Zhao. Orthogonal matrix retrieval with spatial consensus for 3D unknownview tomography. arXiv preprint arXiv:2207.02985, 2022.
 [12] A. Iudin, P. K. Korir, J. SalavertTorres, G. J. Kleywegt, and A. Patwardhan. EMPIAR: a public archive for raw electron microscopy image data. Nature methods, 13(5):387–388, 2016.
 [13] Z. Kam. The reconstruction of structure from electron micrographs of randomly oriented particles. Journal of Theoretical Biology, 82(1):15–39, 1980.
 [14] B. Landa and Y. Shkolnisky. Approximation scheme for essentially bandlimited and spaceconcentrated functions on a disk. Applied and Computational Harmonic Analysis, 43(3):381–403, 2017.
 [15] B. Landa and Y. Shkolnisky. Steerable principal components for spacefrequency localized images. SIAM Journal on Imaging Sciences, 10(2):508–534, 2017.
 [16] C. L. Lawson, A. Patwardhan, M. L. Baker, C. Hryc, E. S. Garcia, B. P. Hudson, I. Lagerstedt, S. J. Ludtke, G. Pintilie, R. Sala, et al. EMDataBank unified data resource for 3DEM. Nucleic Acids Research, 44(D1):D396–D403, 2016.
 [17] C.H. Lee and R. MacKinnon. Structures of the human HCN1 hyperpolarizationactivated channel. Cell, 168(12):111–120, 2017.
 [18] J.Y. Lee and L. Greengard. The type 3 nonuniform FFT and its applications. Journal of Computational Physics, 206(1):1–5, 2005.
 [19] E. Levin, T. Bendory, N. Boumal, J. Kileel, and A. Singer. 3D ab initio modeling in cryoEM by autocorrelation analysis. In 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018), pages 1569–1573. IEEE, 2018.

[20]
W. Liu and J. Frank.
Estimation of variance distribution in threedimensional reconstruction. i. theory.
JOSA A, 12(12):2615–2627, 1995.  [21] N. F. Marshall, O. Mickelin, and A. Singer. Fast expansion into harmonics on the disk: a steerable basis with fast radial convolutions. arXiv preprint arXiv:2207.13674, 2022.
 [22] P. A. Penczek, M. Kimmel, and C. M. Spahn. Identifying conformational states of macromolecules by eigenanalysis of resampled cryoem images. Structure, 19(11):1582–1590, 2011.
 [23] A. Rangan, M. Spivak, J. Andén, and A. Barnett. Factorization of the translation kernel for fast rigid image alignment. Inverse Problems, 36(2):024001, 2020.
 [24] N. Sharon, J. Kileel, Y. Khoo, B. Landa, and A. Singer. Method of moments for 3D single particle ab initio modeling with nonuniform distribution of viewing angles. Inverse Problems, 36(4):044003, 2020.
 [25] Y. Shi and A. Singer. Abinitio Contrast Estimation and Denoising of CryoEM Images. Computer Methods and Programs in Biomedicine, 224:107018, 2022.

[26]
M. Van Heel.
Multivariate statistical classification of noisy images (randomly oriented biological macromolecules).
Ultramicroscopy, 13(12):165–183, 1984.  [27] M. Van Heel and J. Frank. Use of multivariates statistics in analysing the images of biological macromolecules. Ultramicroscopy, 6(1):187–194, 1981.
 [28] M. van Heel, R. Portugal, and M. Schatz. Multivariate statistical analysis in single particle (Cryo) electron microscopy. An electronic textbook: electron microscopy in life science. 3DEM Network of Excellence, London, 2009.
 [29] M. van Heel, R. V. Portugal, M. Schatz, et al. Multivariate statistical analysis of large datasets: single particle electron microscopy. Open Journal of Statistics, 6(04):701, 2016.
 [30] W. Wong, X.c. Bai, A. Brown, I. S. Fernandez, E. Hanssen, M. Condron, Y. H. Tan, J. Baum, and S. H. Scheres. CryoEM structure of the Plasmodium falciparum 80S ribosome bound to the antiprotozoan drug emetine. Elife, 3:e03080, 2014.
 [31] Z. Zhao, Y. Shkolnisky, and A. Singer. Fast steerable principal component analysis. IEEE Transactions on Computational Imaging, 2(1):1–12, 2016.
 [32] Z. Zhao and A. Singer. Rotationally invariant image representation for viewing direction classification in cryoEM. Journal of structural biology, 186(1):153–166, 2014.