The symmetries of image formation by scattering. II. Applications

We show that the symmetries of image formation by scattering enable graph-theoretic manifold-embedding techniques to extract structural and timing information from simulated and experimental snapshots at extremely low signal. The approach constitutes a physically-based, computationally efficient, and noise-robust route to analyzing the large and varied datasets generated by existing and emerging methods for studying structure and dynamics by scattering. We demonstrate three-dimensional structure recovery from X-ray diffraction and cryo-electron microscope image snapshots of unknown orientation, the latter at 12 times lower dose than currently in use. We also show that ultra-low-signal, random sightings of dynamically evolving systems can be sequenced into high quality movies to reveal their evolution. Our approach offers a route to recovering timing information in time-resolved experiments, and extracting 3D movies from two-dimensional random sightings of dynamic systems.


I. INTRODUCTION
In an earlier paper [1], hereafter referred to as Paper I, we presented a theoretical framework for analyzing snapshots formed by scattering. In this paper, we demonstrate the power of this approach to reconstruct threedimensional (3D) models and time-series from random sightings at extremely low signal, with no orientational or timing information. The theoretical framework in Paper I represents the information content of an ensemble of snapshots as a Riemannian manifold, and shows that the properties of operations in space give rise to object-independent symmetries. Purposeful navigation on this manifold is tantamount to reconstructing a 3D model of the sighted system and/or its evolution, in the sense that given any snapshot, any other can be produced on demand. The symmetries of the manifold reveal its natural eigenfunctions, thus allowing physicallybased interpretation of graph-theoretic analysis, and enhanced noise discrimination. Simple algorithms then suffice to reach exceptionally low signal-to-noise levels unmatched by other approaches in terms of computational cost, noise robustness, or both. As examples, we demonstrate structure recovery from radiation-sensitive objects at doses at least an order of magnitude below current levels (signal-to-noise ratio (SNR): −16 dB), and reconstruction of time-series at SNR values as low as −21 dB. The versatility of the approach is demonstrated in the context of simulated and experimental data from X-ray diffraction, cryo-electron microscopy (cryo-EM), and optical snapshots using a variety of graph-theoretic techniques. These applications demonstrate the generality of the symmetry-based approach, elucidating at the same time the measures needed to deal successfully with experimental data, a key benchmark of the practical utility of any theoretical framework. This paper is organized as follows. Without claim to be comprehensive, Sec. II briefly summarizes previous work in the field to provide a context for the applications discussed in this paper. For the convenience of non-mathematical readers, Sec. III provides a conceptual outline of the theoretical framework developed in Paper I. Sec. IV A describes 3D reconstruction from simulated diffraction snapshots of single biomolecules at the signal level expected from single molecules in upcoming experiments utilizing the new generation of X-ray Free Electron Lasers (XFELs) [2][3][4]. Sec. IV B establishes, in principle, the applicability of our approach to crystalline samples. Sec. IV C addresses structure recovery from simulated and experimental cryo-EM snapshots of single molecules. In this case, essential experimental issues such as defocus variation must be faced and incorporated into the theoretical formalism. Sec. IV D demonstrates reconstruction of time-series (movies) from random sequences of ultralow-signal optical snapshots. The paper concludes in Sec. V with a summary of our key findings and their implications. Detailed points of a technical nature are elucidated in appendices, and movies provided as supplementary online material EPAPS [5].

II. PREVIOUS WORK
As described in Paper I, we are concerned with constructing a model from random sightings of a system viewed in some projection, i.e., by accessing a limited number of variables describing the state of the system. A 3D model of an object and its evolution, for example, can be constructed from an ensemble of low-signal 2D snapshots without orientational information [3,4,[6][7][8][9]. Modern graph-theoretic algorithms can now be used to discover low-dimensional manifolds representing the information content of datasets in some high-dimensional space determined by the measurement apparatus [10][11][12][13][14][15][16][17]. The power of these methods stems from their generality, in the sense that no assumptions are made as to the nature of the data or their internal correlations. This brings with it four major challenges: (1) Interpretation of the analysis results ("what physical variables do the manifold dimensions represent?"); (2) Computational cost and scaling behavior on moving from simulated ("toy") datasets to experimental measurements; (3) Robustness against noise, particularly of non-additive, non-Gaussian types; and (4) Incorporation of inevitable and/or desirable experimental factors ("utility of the theoretical framework in practice".) These issues can be brought to focus in the context of the much-discussed problem of recovering the orientation of cryo-EM snapshots of faint biological objects. Direct graph-theoretic attempts to determine the orientation of snapshots from a synthetic object were abandoned at a signal-to-noise ratio (SNR) of ∼ 2 dB, even though only additive Gaussian noise was included [18]. Noting that graph-theoretic analyses often "fail to solve the cryo-EM problem, because the reduced coordinate system that each of them obtains does not agree with the projection directions" [16], properties specific to cryo-EM images were used to extract information from the snapshots. Graphs were then constructed using this information in order to assign physical meaning to the outcome of the analysis. Orienting low-signal cryo-EM snapshots by utilizing so-called (straight) common-lines identified primarily in simulated data with additive Gaussian noise has reached remarkably low SNR values [19]. However, such assumptions, while justified under some circumstances, are not generally valid. Common-lines, for example, are present only when elastic single-scattering dominates, are straight only when the wavelength of the incident radiation is so short that the Ewald sphere can be replaced by a plane, and are compromised by defocus variations essential for reliable structure recovery by cryo-EM [20].
Symmetry-based assignment of physical meaning to the outcome of graph-theoretic analysis of scattering data and its favorable computational consequences were addressed in Paper I. Here, we are concerned with ability of this theoretical framework to deal with noise and other important factors encountered in experimental datasets. This determines the practical utility of an approach as much as theoretical elegance and computa-tional efficiency. Below, we demonstrate the utility of our symmetry-based approach by applying a number of manifold-embedding techniques to a variety of simulated and experimental datasets (see Table I).

III. CONCEPTUAL SUMMARY OF THEORETICAL FRAMEWORK
A snapshot formed on a 2D detector by scattered radiation from an object can be represented by a vector, with the intensity values recorded at the n detector pixels as components (Paper I, Fig. 1). Object motion and/or evolution (dynamics) change the pixel intensities, causing "the vector tips" representing the ensemble of snapshots to trace out a surface -a manifold -in the ndimensional data space. The number of degrees of freedom available to the object determines the dimensionality of the manifold traced out. Rotations of a rigid object in 3D, for example, result in a 3D manifold.
The data manifold represents the totality of information about the object gathered by the detector in the course of an observation. Learning is tantamount to understanding the properties of this manifold sufficiently to "navigate" on it. Learning the manifold generated by object rotations, for example, is equivalent to constructing a 3D model of the object, because, starting from any 2D projection (point on the manifold) any other 2D projection can be found by navigation, with the shortest route corresponding to a geodesic.
As we are initially concerned with constructing 3D models from 2D snapshots, we consider a formulation of scattering by a single object in Fourier space so as to concentrate on the effect of rotations. This circumvents issues such as rigid shifts, which would otherwise have to be corrected or incorporated as additional manifold dimensions. Rotation operations do not commute. One must therefore consider the order in which they are performed. This leads to a distinction between so-called left and right "translations," where a rotation operator T is placed to the left or right of another rotation operator R, i.e., TR vs. RT. A left translation can be thought of as an active rotation in 3D space of the incident beam-detector arrangement (frame rotation) after R. Similarly, a right translation corresponds to an active rotation of the object (object rotation) before R. As TR = RT, left and right translations must be considered separately. Each forms an SO(3) group, and the total set of possible operations to be considered corresponds to SO(3) × SO (3).
A key question is this: Which, if any, of these operations leave the distances on the manifold unchanged, i.e., which operations are "invisible" to an ant crawling on the manifold? These operations would represent symmetries -more precisely isometries -of the manifold. For a detector with circular symmetry, the distances on the manifold are invariant under beam-detector rotations about the beam axis. This is obvious; a frame rotation about the beam axis rotates all the snapshots by the same amount about that axis without changing them. This leaves the distances on the manifold unchanged. The process of image formation on a circularly-symmetric detector at right angles to the illuminating beam thus has SO(2) isometry, i.e., of all possible SO(3) frame operations, the SO(2) subset of rotations about the beam direction leave the distances on the manifold unchanged. This is related to the projection of a 3D object on the 2D detector, which is equivalent to a "central slice" through the diffraction volume in reciprocal space.
Consider next the SO(3) set of operations corresponding to object rotations. It turns out that the metric measuring distance on the manifold can be decomposed into a homogeneous part, which varies uniformly with object rotation, plus a residual term, which acts as a fingerprint of the object (see Paper I Sec. III C). Considering the homogeneous part only, the total set of symmetries, is then SO(2) × SO(3). The same set of symmetries appears in certain models of the universe in general relativity [21,22], and is associated with well-known eigenfunctions familiar in the context of spinning tops in classical and quantum mechanics [23].
The key point here is that the knowledge of the manifold symmetries, which stem from the nature of operations in space, allows one to determine the leading-order properties of the manifold under a very general set of scattering scenarios, including its natural eigenfunctions. Projection of noisy datasets on these eigenfunctions is tantamount to noise discrimination. The components of a data point representing a snapshot can then be directly related to its orientation (see Paper I Sec. III D).

IV. APPLICATIONS
It has long been known that the use of problem-specific constraints can substantially increase computational efficiency [24]. By combining wide applicability with class specificity, symmetries represent a particularly powerful example of such constraints. In Paper I, we used the object-independent symmetries of image formation to recover 3D structure from a large ensemble of simulated, noise-free diffraction snapshots with a computational complexity 10 4 × higher than the state of the art. Here we demonstrate the noise robustness stemming from exploiting the symmetries of image formation. Examples include orientation recovery, 3D reconstruction, and movie extraction from ultra-low-signal diffraction or image snapshots of periodic and non-periodic objects and dynamical systems. Each example was selected to highlight an important application area. As shown in Table I, a variety of manifold-embedding techniques can be used.
A. Structure recovery from simulated diffraction snapshots of non-periodic objects at ultra-low signal First, we demonstrate 3D structure recovery from a collection of two million simulated diffraction snapshots of the synthetic protein chignolin (Protein Data Bank (PDB) descriptor: 1UAO, model 1) at 4 × 10 −2 scattered photons per Shannon pixel at 1.8Å. (A Shannon pixel is of the size needed for appropriate sampling of the intensity distribution as prescribed by the Shannon-Nyquist theorem.) This scattered intensity is expected from a 500 kD protein exposed to a single pulse from an XFEL [2,3]. At this signal level, Poisson (shot) noise dominates. The ability to deal with such levels of non-additive noise was previously demonstrated only by Bayesian algorithms [3,8] with extremely unfavorable scaling behavior [1,3,4,8], restricting the size of amenable objects to eight times the spatial resolution.
Here, we use the symmetry-based approach described in Paper I after modest denoising. The denoising scheme consists of two steps: (1) Convolve the snapshot pixels with a 2D Gaussian filter with a width approximately equal to that of a Shannon pixel; (2) Replace each snapshot vector by an average over its local neighbors. Depending on the SNR, a number iterations of step (2) may be needed, with a stopping criterion based on a leastsquares residual determined through the first nine Laplacian eigenfunctions of the dataset (ordered in order of increasing eigenvalue). These eigenfunctions are employed in our scheme to assign an orientation to each snapshot. For details see Appendix A and Paper I.
To estimate the accuracy of orientation recovery, we use the following measure for root-mean-square (RMS) distance between the deduced and true orientations: where D ij = 2 arccos(|τ i · τ j |) andD ij = 2 arccos(|τ i ·τ j |) are the true and estimated internal distances between orientations i and j, respectively, and · is the inner product between quaternions. Moreover, to assess the influence of local averaging on the eigenfunctions employed for orientation recovery, we compute the distance γ of the invariant subspaceṼ spanned by the leading nine eigenfunctions of the diffusion matrix P in Table V from the corresponding invariant subspace V of the noise-free diffusion matrix [25]. Note that P has size s × s, where s is the number of snapshots in the data set; i.e.,Ṽ and V are subspaces of R s . Here, we employ a standard distance measure from matrix perturbation theory [26], viz.
whereΠ and Π are orthogonal projectors from R s toṼ and V , respectively, and · 2 denotes the spectral norm of matrices. With this definition, γ lies in the interval [0, 1], and may be interpreted as the sine of an angle characterizing the deviation ofṼ from V . For our purposes, Eq. (2) is more appropriate than an error measure based on the difference between the noisy and noise-free diffusion matrices (or their generators), since the latter depends on higher eigenfunctions which are not used in our scheme. Diffraction snapshots were simulated in 2 × 10 6 different orientations to a spatial resolution of 1.8Å using 1Å photons. The orientations were sampled approximately uniformly over SO(3), as described in Ref. [27]. Cromer-Mann atomic scattering factors [28] were used for the 77 non-hydrogen atoms, and the hydrogen atoms neglected. The detector pixel was the appropriate Shannon pixel [3], which oversamples the scattered amplitudes by a factor of two, resulting in 40 × 40 = 1600 Shannon detector pixels. To model shot noise, diffracted intensities were scaled so that the mean photon count (MPC) per Shannon pixel was 0.04 at 1.8Å resolution. The quantized photon count at each pixel was obtained from a Poisson distribution by the algorithm described in Ref. [29].
With no other information, the noisy diffraction patterns were provided to the algorithm in Table III (width of Gaussian filter for image smoothing σ = 0.7; number of nearest neighbors in the sparse distance matrix d = 220; number of nearest neighbors for local averaging l = 20; number of datapoints for least-squares fitting r = 8 × 10 4 ; number of nearest neighbors for autotuning n = 30.) As illustrated in Fig. 1, the least-squares residual G * , the subspace distance γ, and the RMS orientation recovery error ε all decrease monotonically for the first five iterations of local averaging, but exhibit a mild increase at iteration six. At that point the algorithm was terminated in accordance with the stopping criterion described above and in Appendix A. The minimum ε value attained with this choice of parameters at iteration 5 is ∼ 1.1 Shannon angles. We measured comparable levels of orientation-recovery accuracy for various combinations of l and n parameters in the range 10-50. In all cases, we observed that small values of G * correlate strongly with small values of , indicating that the least-squares residual provides an effective guideline for setting the parameters of the algorithm. This is particularly important, because G * depends solely on the Laplacian eigenfunctions (see Table IV), and, unlike ε, can be evaluated in an experimental environment where the correct orientations are not known.
The quality of orientation recovery was further tested by inverting the reconstructed 3D diffraction volume compiled on a uniform Cartesian grid by an interpolation scheme consistent with the geometry of diffraction [30]. The R-factor between the gridded scattering amplitudes F i and those obtained from the Fourier transform of the recovered electron density from phasing, F i , was defined as The 3D electron density obtained by iterative phasing with the superflip algorithm [31] (R = 0.20) is shown in Fig. 2. The close agreement with the known structure of chignolin clearly demonstrates sufficient alignment accuracy for reconstruction to 1.8Å resolution. This is on par with the computationally much more expensive Bayesian approaches [3,8,32].

B. Orienting diffraction patterns of crystals
So far we have shown that our symmetry-based approach can be used to orient diffraction patterns from single molecules to high accuracy. We now demonstrate that this approach can also orient diffraction snapshots from crystals. This is important, because recent XFEL-based "diffract-and-destroy" approaches, which use femtosecond X-ray pulses to "outrun radiation damage", produce diffraction snapshots of nanocrystals of unknown orientation [33]. As a representative example, we consider a biological crystal of the enzyme superoxide dismutase-1 (SOD1, PDB designation: 1AR4) with ∼ 3 × 10 3 atoms per unit cell, and thus highly complex diffraction patterns. The key issue is whether manifolds produced by diffraction snapshots of crystals are sufficiently homogeneous (possess sufficiently homogeneous metrics) for snapshot orientations to be recovered in a straightforward manner. To demonstrate this point, we intentionally utilized snapshots spanning an orientation range of 90 • so as to produce an open 1D manifold, and analyzed the dataset with the Isomap manifold-embedding method [10]. In contrast to Diffusion Map, whose eigenfunctions are insensitive at the boundaries, Isomap maps a 1D open manifold to a straight line segment, and is sensitive to snapshot orientation over the entire range.
Experimental diffraction patterns of a single crystal of superoxide dismutase-1 with a mosaicity of 0.8 • were obtained at the Advanced Photon Source (λ = 0.98Å). The crystal was rotated about an arbitrary axis with a step size of 0.5 • , and 1800 diffraction patterns recorded over a range of 90 • . To compensate for spurious beam intensity fluctuations, the diffraction pattern intensities were normalized. Isomap was used to embed diffracted amplitudes (square-roots of intensities), using two nearest neighbors for calculation of geodesic distances (integrals of the metric). As shown in Fig. 3, a onedimensional and uniformly populated manifold results, with the projection on the first eigenvector linearly proportional to the snapshot orientation to within 1 • , compared with the crystal mosaicity of 0.8 • . The homogeneity of this manifold (metric) establishes that, in principle, our symmetry-based approach can be used to treat crystalline objects in the same way as non-periodic single particles, provided, of course, object symmetry is appropriately incorporated.

C. Structure recovery from experimental cryo-electron micrographs
A well-studied application of graph-theoretic techniques concerns the 3D reconstruction of faint biological objects by single-particle cryo-EM without orientational information. In cryo-EM, the resolution is strongly degraded by radiation damage. As such, the lowest acceptable exposure to electrons and thus SNR must be employed. As described in Sec. II, this has proved a fertile ground for testing new algorithms. By recourse to specific properties of cryo-EM images, impressive results have been obtained, primarily with simulated data. Beyond noise, however, reconstruction by cryo-EM must contend with a range of key issues, chief among them the loss of information due to zero-crossings in the transfer function of the microscope and thus partial loss of information in any single snapshot. The exact position of the zero-crossings depends sensitively on microscope defocus. This offers a means to recoup some of the lost informa-tion by insuring that the dataset includes micrographs obtained over a range of defocus values, each with a different set of zero-crossings in the transfer function. The key point is this: the object structure cannot be recovered in full detail from a single defocus, even if the imaging parameters were known exactly. Thus, for a reconstruction algorithm to be of practical use, it must deal with the effect of defocus variations -a test rarely passed by new algorithms. Here, we demonstrate structure recovery from experimental cryo-EM images of the biological molecule chaperonin. Specifically, we incorporate the effect of defocus, use the symmetry-based homogeneity of the manifold metric to deduce orientations, and thence recover the 3D object structure. This demonstration is mitigated by two factors: (1) in order to expedite the calculations, snapshots with only one orientational degree of freedom were selected from a set presorted by a standard orientation algorithm; and (2) to demonstrate structure recovery at ultra-low signal -far below what is normally used -experimental snapshots were preprocessed to simulate such low signal levels.
Randomly oriented single-particle cryo-EM images images of the wild-type group II chaperonin in methanococcus maripaludis (Mm-cpn), obtained with a mean incident electron count (MEC) of 20/Å 2 (equivalent to 135 electrons per 2.6Å-square snapshot pixel) were kindly provided by Chiu et al. [34]. Each snapshot consisted of 96 × 96 pixels. A set of 413 side-view snapshots was selected from 5000 images, whose orientations had been previously determined by the eman program [35], resulting in a dataset with a single orientational degree of freedom about the object symmetry axis.
To investigate the performance of our method at lower dose, a second data set was produced by applying an additional Poisson process to the raw experimental images. The method is based on an approximation valid for lowcontrast images with Poisson noise and sufficiently large MEC. The substitution I → I = Pois(I 1/2 ) transforms a signal I to a signal I , with mean MEC = MEC 1/2 and variance var(I ) = var(I)/4. Simulations verified the accurate validity of this approach at MEC = 100, compared with an MEC per snapshot pixel of 135 for our experimental images. Twenty noisy versions of each image were thus generated to form a data set of 8260 images with an effective MEC of 1.7 perÅ 2 .
Since neither the noise-free signal nor the noise variance was known for our experimental cryo-EM images, a method developed by Frank [20] was used to estimate the SNR directly from the experimental data. This determines the SNR from the cross-correlation coefficient C ij between two images in the same orientational class using the definition: SNR = 10 log 10 mean(C ij /(1 − C ij )), where the mean is taken over all classes and all images within each class. Provided two images represent different realizations of noise from an identical object in the same orientation, the above estimate for SNR agrees with the standard definition SNR = 10 log 10 var(signal) var(noise) [36]. With the classification obtained from eman [35], and the assumption that class members differ only in noise, we estimate a SNR of −6 dB for the raw experimental snapshots (MEC: 20/Å 2 ) and a SNR of −16 dB for the preprocessed experimental snapshots (MEC: 1.7/Å 2 ).
As described above, the defocus value and hence transfer function of the electron microscope vary from snapshot to snapshot. To analyze such cryo-EM data, we implemented a modified version of the manifoldembedding algorithm Generative Topographic Mapping (GTM) [17,37] to explicitly incorporate the effect of the microscope transfer function. GTM defines a manifold in data space by partitioning the noisy dataset into a number of Gaussians each centered around a point (node) on the manifold. The partitioning is based on a nonlinear mapping of a latent space, in this case the space of rotations. GTM is thus, in essence, a manifold-embedding technique, with the symmetries of scattering manifested in the homogeneity of the data manifold, as described in Paper I. However, the generative capability of GTM allows one to construct an image (in essence a model snapshot) at each node on the data manifold. In our approach, this model image extracted from the data corresponds to the aberration-free projected potential of the object. In order to assign an experimental snapshot to a model image, its distance from the model is calculated after convolving the model with the transfer function of the microscope at the defocus corresponding to that of the experimental snapshot. This convolution proceeds efficiently as multiplication in Fourier space, and is not computationally expensive. A similar approach based on more efficient manifold-embedding techniques will be published elsewhere.
The GTM-based approached was first validated with simulated cryo-EM images of chaperonin over a typical experimental defocus range of 10, 000Å to 30, 000Å (underfocus). The orientations were successfully recovered to within 1 • . To reconstruct 3D density maps from experimental snapshots, model aberration-free projected potentials were generated (lifted) at 16 equally-spaced nodes of the data manifold produced by experimental images replicated according to the 8-fold object symmetry. 3D density maps were then reconstructed tomographically using the back-projection algorithm bg cg of the spider software package [38]. For comparison, a simulated density map was obtained from 2D snapshots using the known chaperonin atomic coordinates (PDB identifier: 3LOS) under the following imaging conditions: spherical aberration C s = 4.1 mm; defocus ∆f = 24,000Å (underfocus); electron energy E = 300 keV; damping envelope parameter B = 50Å 2 ; images phase-flipped. The resulting 3D density map was passed through a 5Å Gaussian filter. Fig. 4(a) shows a typical experimental snapshot, Fig. 4(b) the average of the micrographs assigned to an orientation class by the cryo-EM reconstruction software package eman [35], and Fig. 4(c) the snapshots oriented by manifold embedding and reconstructed (lifted) from the manifold. (For a movie of the reconstructed tilt series see EPAPS Movie 1). Note that the manifold is able to generate missing images by interpolation. The improved quality of the manifold-generated snapshots compared to the class averages offers the possibility to reconstruct at significantly reduced dose. Fig. 4(d) is an experimental snapshot preprocessed to approximate snapshots expected from a single chaperonin molecule at a dose 12× lower than commonly used [34] (SNR ∼ −16 dB, i.e., 10 dB below a typical dose). Fig. 4(e) is the snapshot lifted from the manifold after orienting an ensemble of 8000 different raw snapshots by manifold embedding. It is clear from this image, the corresponding tilt series (EPAPS Movie 2), and the 3D reconstructions of Fig. 4(f-h) that snapshots can be successfully oriented by manifold embedding to produce 3D models, even at 12× lower signal than in use today. Note that images at this dose could not be oriented by standard cryo-EM approaches [38], even when accurately centered prior to analysis, as was performed here. In contrast, our orientation recovery results were similar to those obtained at an MEC of 20/Å 2 , indicating that the effect of lower signal levels can be compensated by increasing the number of snapshots. Our results thus offer the tantalizing possibility of reducing the snapshot dose in 3D reconstruction techniques using ionizing radiation, in some cases by at least an order of magnitude. This would significantly mitigate the limits set by radiation damage. As a benchmark, the essentially unfulfilled promise of the costly transition of cryo-EM to liquid He temperatures was aimed at improving dose tolerance by a factor of two. The superior signal extraction capability offered by manifold mapping could also be used to obtain images at smaller defocus values in order to reconstruct the object to higher resolution.

D. Time-series (movies) from ultra-low-signal random-sequence snapshots
Our knowledge of the precise time at which an experimental snapshot of a dynamic system was obtained is corrupted by inevitable uncertainties, which can substantially exceed the intrinsic time resolution of the observation technique. Modern pump-probe experiments, for example, can now be performed with pulses as short as a few femtoseconds, but their time-resolution is often determined by timing jitter, which can be up to two orders of magnitude larger [39,40]. When the state of a system under observation is not synchronized with the observation windows, a sequence of snapshots can represent random sightings of the system during its evolution. It is thus important to develop means for deducing "time stamps" directly from the data, either to reduce jitterinduced uncertainty, or to order a sequence of snapshots according to the intrinsic evolution of the system un- der observation. Here, we demonstrate this capability at SNRs as low as −21 dB. Specifically, we show that: (1) Randomized movie sequences can be time-ordered, even when the signal is extremely low; and (2) Frames generated (lifted) from the manifold produce movies of superior quality. Movies of a pirouette and a pas de deux were downloaded from the web. These represent optical snapshots of a conformationally rigid body in rotations, and the evolution of two flexible bodies in interaction, respectively. In order to reduce the SNR, a constant background was added, and shot noise incorporated at each pixel depending on its intensity value, as described in Ref. [29]. For the pirouette, a sequence of 16 turns consisting of 268 frames (210×160 pixels each) was replicated 132 times, a background 5× the mean intensity added, and shot noise incorporated to produce an effective mean photon count per pixel of 0.08 and a SNR of −21 dB (see Eq. 4). For the pas de deux, a sequence of 870 frames (265 × 305 pixels each) was replicated 12 times with an added background of twice the mean intensity, and shot  4: Initialize a set of weights W using principal component analysis. 5: Initialize inverse Gaussian noise variances α and β. 6: Compute a set of responsibilities R that assigns each snapshot to a node from the results of manifold embedding.
8: repeat 9: where the regularization parameter λ may be zero.
13: until convergence of W 14: M ← ΦW 15: return M noise incorporated to produce a mean photon count of 0.8 and a SNR of −11 dB. For both movies, camera motion was corrected by reference to a stationary marker.
Each random sequence was ordered by a suitable manifold-embedding technique (Diffusion Map or Isomap). Using the generative property of GTM, images were then lifted from the manifold. As described in Ref. [4] and demonstrated in Fig. 4(c,e), this procedure uses the information content of the entire dataset to generate each snapshot, producing images of significantly higher quality than possible by traditional classifying and averaging techniques. It is also more robust against nonuniform sampling and jitter. Table II summarizes the  lifting procedure. For the pirouette, Diffusion Map was used to recover the object orientation in each frame during the dance, and hence the sequence order (number of nearest neighbors in the sparse distance matrix d = 5896; Gaussian kernel bandwidth = 200.) Snapshots were then lifted from the manifold (number of nodes K = 28; number of basis functions M = 14, basis function width σ = 2.) Reconstructed images are shown in Fig. 5 together with a sequence of unsorted, unprocessed snapshots. The movie is available as EPAPS Movie 3. The randomized pas de deux sequence was ordered with Isomap (number of nearest neighbors d = 33.) To compile the movie, images were lifted from an ordered sequence of 870 points (nodes), corresponding to uniform sampling on the Isomap manifold. Reconstructed images are shown in Fig. 6 together with a sequence of unsorted snapshots. The movie is available as EPAPS Movie 4.
Figs. 5 and 6, and the associated movies clearly show that our approach is able to determine the correct frame sequence and generate high quality snapshots at signal levels as low as 0.08 photon/pixel for the pirouette (modulo one revolution), and 0.8 photon/pixel for the pas de deux, both with added background and non-additive noise corresponding to signal-to-noise ratios in the range −11 to −21 dB. These examples demonstrate the capability to determine the time evolution of systems from unsorted random sightings at extremely low signal. They also highlight the potential to correct timing jitter in pump-probe experiments, and reconstruct the evolution of dynamic systems from random sightings of members of a heterogeneous set, each at a different stage of its evolution. These possibilities will be described in detail elsewhere. The general implications for signal extraction and image processing are clear.

V. CONCLUSIONS
We have shown that manifold mapping, as described in Paper I, augmented with modest noise reduction measures, is able to extract structural and timing information from simulated and experimental snapshots at extremely low signal. The ability to orient simulated and experimental diffraction and image snapshots confirms the accessibility of the homogeneous manifold expected from our theoretical framework for a wide range of objects and imaging scenarios, including crystalline samples. The capability to recover 3D structure at extremely low signal is on par with the more expensive Bayesian approaches, but offers greater reach in terms of sample size and resolution, as demonstrated in Paper I. The noiserobustness of our approach substantially exceeds what has been demonstrated with comparable graph-theoretic approaches without restrictive, application-specific assumptions. The manifold itself offers a powerful route to image reconstruction at low signal, because snapshots reconstructed from the manifold achieve higher signal-tonoise ratios than possible by traditional approaches based on classification and averaging. Taken together, these offer a physically-based, computationally efficient, noiserobust route to analyzing the large and varied datasets generated by existing and emerging structure recovery methods. In the longer term, it should be possible to use these approaches to recover or improve timing information in pump-probe experiments, and construct 3D movies (4D maps) from random sightings of members of structurally heterogeneous and dynamically evolving ensembles.

ACKNOWLEDGMENTS
We are grateful for experimental data and advice to W. Chiu and J. Zhang (cryo-EM images) and M. Schmidt (crystallographic data), and acknowledge discussions with G. N. Phillips, Jr., R. Rosner, W. Schröter, and members of the UWM Physics Department. We are grateful to J. Frank  In the manifold picture, noise can be described as a perturbation of the noise-free manifold M = {a 1 , a 2 , . . . , a s }, where a i is a snapshot vector of mea-sured pixel amplitudes, viz.
This causes the observed, noisy, dataset not to lie exactly on the manifold M (up to a global scaling by κ). One would expect that if the perturbation norm δa i becomes comparable to the kernel bandwidth 1/2 , the computed eigenfunctions ψ k are distorted to the point that the embedded manifold no longer has the topology of SO(3) [15,41]. Indeed, as illustrated in Fig. 1 direct application of the algorithm for noisefree data (see Paper I Table III) to a noisy datasetM at an MPC = O(10 −2 ) results in poor accuracy. In order to be practically useful, the noise-free orientationrecovery scheme must be augmented by a suitable denoising method. Conceptually, we denoise the data in three steps: (1) Low-pass filtering of each snapshot by convolution of pixel intensities with a 2D Gaussian; (2) Variancestabilizing transformation (VST); and (3) Local averaging over nearest neighbors in data space prior to embedding. In practice, we combine (1) and (2), known to be effective for shot noise [42][43][44], into a single step, and follow the iterative procedure described below.

Low-pass Gaussian filtering and variance-stabilizing transformation (VST)
Convolution with a low-pass Gaussian filter of bandwidth σ is represented by: with VST, proposed by Guan [44] for low-intensity data, can be written as: Noise robustness can be further enhanced by a so-called self-tuning Gaussian kernel introduced by Zelnik-Manor and Perona [45]. Here, instead of the isotropic Gaussian kernel K(a i , a j ) = exp(− a i − a j 2 / ) of Eq. (B3) in Paper I, one uses an anisotropic Gaussian kernel with local scaling parameters, i , given by A canonical choice for the scaling parameters, which we adopt throughout, is i = a i − a N il 2 , where, as usual, N il denotes the index of the l-th nearest neighbor of datapoint i.

Iterative local averaging
If the true nearest neighbors of a point in data space are known, and noise produces no systematic bias, local averaging approaches the true manifold. To see this, let ε be an error tolerance in data space, and consider a reference orientation R with corresponding noise-free snapshot a. For any ε > 0 it is possible to find a ball B ε in data space centered at a, such that for any countable set {a 1 , . . . , a l } of noise-free snapshots lying in B ε the mean,ā = l i=1 a i /l, has error ā − a < ε. In the presence of noise, the a i are replaced by the random variables in Eq. (A1a); i.e., a i →ã i , whereã i are statistically independent, have expectation value κ 1/2 a i proportional to a i , and finite variance ∆a 2 i . Moreover, the sample mean within B ε becomes a random variableâ = l i=1ãi /l with expectation value κ 1/2ā and variance ∆a 2 = l i=1 ∆a 2 i /l. By the law of large numbers, in the limit of an infinite data set with infinite snapshots in B ε (i.e., l → ∞),â i is equal toā i (up to an unimportant proportionality constant) with probability one. Thus, recovery of the data manifold with error ε is possible almost surely.
In practice, noise corrupts the local neighborhood relations, and without a priori information, it is not possible to identify which of the snapshots in a noisy data set are associated with the ball B ε of the underlying noise-free system. Therefore, it cannot be guaranteed that local averaging leads to the correct manifold. We therefore exploit our knowledge of the natural eigenfunctions of scattering manifolds to monitor the effect of local averaging, and terminate the procedure before substantial deviations have occurred.
Specifically, we follow the algorithm described in Table III. First, we apply the VST operation (A5) to the intensity data {I i } = M I , setting the filter width σ to a relatively small value (e.g., in Sec. IV A, σ is set to

10:
Execute the algorithm in Table IV with input dataMi; store the outputs as Ti, Ni, and G * i .

11
: terminate if the residual has increased.

17:
7/10 of a pixel width). The autotuning version of the orientation-recovery method (Table IV) is executed using the VST-filtered intensities as input data. The nearestneighbor indices N 0 obtained in the course of the calculation of the sparse distance matrix then become our initial estimate for the true nearest-neighbor indices. We also record the residual of the nonlinear least-squares output, G * 0 , and choose a value l ≤ d for the number of nearest neighbors for local averaging.
Next, we enter an iteration loop, where in the i-th step the datasetM is computed, and the algorithm in Table IV executed usingM i as input data. The resulting quaternion estimates, nearest-neighbor indices, and least-squares residual are respectively designated T i , N i , and G * i . Note that the residual G * i is a measure of the difference between the eigenfunctions obtained by embedding and the natural eigenfunctions (Wigner D-functions) expected on the basis of symmetry (see Paper I).
If, after an iteration i, G * i is larger than the residual G * i−1 encountered in the previous step, the loop is terminated. Otherwise, the iteration is repeated using N i as an updated estimate of the true nearest-neighbor indices. Our final orientation (quaternion) assignment is the one corresponding to the minimum least-squares 2: return N 3: Rescale the distance data by the n-th nearest neighbors: Sij ← S i,N i,n S j,N j,n 1/2 .

4:
Compute the sparse transition probability matrix P using the algorithm in Table V with inputs S, N, , and α = 1.

10:
ProjectRi to an orthogonal matrix

11:
Convert Ri to a unit quaternion τi 12:   The empirical evidence in Sec. IV A clearly shows that, given a sufficiently large number of sample points, the scheme, applied only a handful of times and terminated using the value of G * i as a criterion, provides noise reduction sufficient for accurate orientation recovery at MPC = O(10 −2 ).
A potentially fruitful way of interpreting mathematically the success of the process (which lies outside the scope of the present paper) would be to explore its connections with mutually reinforcing models (MRMs) for graph filtering [46]. This type of model involves iteratively replacing vertices of graphs with weighted averages, whereby the vertex itself and its local neighborhood exert an influence on the vertex in the course of iterative updates. In certain applications, the iterative process in an MRM is terminated after only a small number of iterations. Both of these two features are present in the scheme presented here.