The Symmetries of Image Formation by Scattering. I. Theoretical Framework

We perceive the world through images formed by scattering. The ability to interpret scattering data mathematically has opened to our scrutiny the constituents of matter, the building blocks of life, and the remotest corners of the universe. Here, we deduce for the first time the fundamental symmetries underlying image formation. Intriguingly, these are similar to those of the anisotropic"Taub universe"' of general relativity, with eigenfunctions closely related to spinning tops in quantum mechanics. This opens the possibility to apply the powerful arsenal of tools developed in two major branches of physics to new problems. We augment these tools with graph-theoretic means to recover the three-dimensional structure of objects from random snapshots of unknown orientation at four orders of magnitude higher complexity than previously demonstrated. Our theoretical framework offers a potential link to recent observations on face perception in higher primates. In a later paper, we demonstrate the recovery of structure and dynamics from ultralow-signal random sightings of systems with no orientational or timing information.


I. INTRODUCTION
We perceive by constructing three-dimensional (3D) models from random sightings of objects in different orientations. Our ability to recognize that we are seeing a profile even when the face is unknown [1] suggests that perception may rely on object-independent properties of the image formation process. The discovery of these properties would underpin our mathematical formalisms, enhance our computational reach, and perhaps elucidate the process of perception.
Here, we show that image formation by scattering possesses specific symmetries, which stem from the nature of operations in three-dimensional (3D) space. This follows from a theoretical framework, which considers the information contained in a collection of random sightings of an object. An example is a collection of twodimensional (2D) snapshots of a moving head, or a rotating molecule. Our primary results can be summarized as follows. (1) The information gleaned from random sightings of an object by scattering onto a 2D detector can be represented as a Riemannian manifold with properties resembling well-known systems in general relativity and quantum mechanics. As shown below and in a subsequent paper, this allows one to efficiently construct a 3D model of the object. (2) The information about an object can be decomposed into an object-independent term reminiscent of a Platonic Form or "ideal object" [2], plus an object-specific fingerprint. (3) The objectindependent term can be used to determine the object * All authors contributed equally. Correspondence should be addressed to: ourmazd@uwm.edu orientation giving rise to each snapshot independently of the object, while the object-specific fingerprint can be used for recognition purposes. This seems to have been experimentally observed in face recognition in higher primates [1].
More generally, it is now well-established that numerical data clouds give rise to manifolds, whose properties can be accessed by powerful graph-theoretic methods [3][4][5][6][7][8]. However, it has proved difficult to assign physical meaning to the results [9,10]. Using the arsenal of tools developed in differential geometry, general relativity, and quantum mechanics, we elucidate the physical meaning of the outcome of graph-theoretic analysis of scattering data, without the need for restrictive a priori assumptions [11]. Finally, perception can be formulated as learning some functions on the observation manifold. Our approach is then tantamount to machine learning with a dictionary acquired from the empirically accessible eigenfunctions of well-known operators.
After a brief conceptual outline in Sec. II, we summarize relevant previous work in Sec. III. The symmetrybased scheme for analyzing scattering data is developed in Sec. IV, and the utility of these concepts demonstrated in Sec. V in the context of simulated images from X-ray scattering. We present our conclusions in Sec. VI. Material of a more technical nature, including mathematical derivations and algorithms, is provided in Appendices A and B. A movie of a reconstructed object in 3D is presented as supplementary online material [12]. Further applications are described in a subsequent paper [13], hereafter referred to as Paper II.

II. CONCEPTUAL OUTLINE
We are concerned with constructing a model from sightings of a system viewed in some projection, i.e., by accessing a subset of the variables describing the state of the system. A 3D model of an object, for example, can be constructed from an ensemble of 2D snapshots. Each snapshot can be represented by a vector with the pixel intensities as components. A collection of snapshots then forms a cloud of points in some high-dimensional data space (Fig. 1). In fact, the cloud defines a hypersurface (manifold) embedded in that space, with its dimensionality determined by the number of degrees of freedom available to the system. Snapshots from a rotating molecule, for example, form a 3D hypersurface.
This perspective naturally leads one to use the tools of differential geometry for data analysis, with the metric playing a particularly important role. In non-technical terms, one would like to relate infinitesimal changes in a given snapshot to the corresponding infinitesimal changes in orientation, giving rise to the changes in the snapshot. In other words, one would like to relate the metric of the data manifold produced by the collection of snapshots to the metric of the manifold of rotations. This would allow one to determine the rotation operation connecting any pair of snapshots. The problem, however, is that the metric of scattering manifolds is not simply related to that of the rotation operations. We solve this problem in two steps. First, we show that the metric of the data manifold can be decomposed into two parts, one with high symmetry, and an object-specific "residual" with low symmetry. Second, using results from general relativity and quantum mechanics, we show that the eigenfunctions of the high-symmetry part are directly related to those of the manifold of rotations. This allows one to deduce the snapshot orientations from the highsymmetry part, which is the same for all objects, and use the object-specific part as a fingerprint of each object.
We conclude this section with four observations. First, real datasets necessarily contain a finite number of observations. As such, they must be treated by graphtheoretic means, which tend to the differential geometric limit under appropriate conditions. This issue is addressed below, as needed. Second, we have not distinguished between image and diffraction snapshots. If the dataset consists of the latter, each snapshot, or a reconstructed 3D diffraction volume must be inverted by socalled phasing algorithms [14][15][16]. As this procedure is well established, we do not address it further. Third, the discussion is restricted to the effect of operations on objects with no symmetry. The case where the object itself has specific symmetries will be treated elsewhere. Finally, the knowledge gained in the course of our analysis is sufficient to navigate from any starting point to any desired destination on the manifold; i.e., given any snapshot, produce any other as required. This is, of course, tantamount to possessing a 3D model of the object. As this approach is somewhat unfamiliar, we provide actual FIG. 1. Intensity values of a snapshot on a detector with n pixels, represented as an n-dimensional vector. Changes in the object or its position alter the snapshot, causing the vector to trace out a manifold in the n-dimensional data space. The dimensionality of manifold is determined by the number of degrees of freedom available to the object. For instance, the rotations of a rigid object give rise to a 3D manifold.
3D models to demonstrate the power of our approach.

III. PREVIOUS WORK
In order to place our contribution in context, we present a brief summary of previous work. As a review is beyond the scope this paper, the summary is necessarily brief, leaving out much excellent work in this general area.
When the snapshots emanate from known object orientations and the signal is adequate, standard tomographic approaches [17] are routinely employed. 3D models can also be constructed when the snapshot orientations are unknown and the signal is low [18] [signal-to-noise ratio (SNR) ∼ −5 dB], or extremely low [19][20][21] (∼ 10 −2 scattered photons per Shannon pixel with Poisson noise and background scattering). Finally, efforts are underway to recover structure, map conformations, and determine dynamics ("3D movies") from random sightings of non-identical members of heterogeneous ensembles [21][22][23] and/or evolving systems [21,23].
Data-analytical approaches now include powerful graph-theoretic and/or differential geometric means to deduce information from the global structure of the data representing the (nonlinear) correlations in the dataset [3-8, 24, 25]. Ideally, this structure takes the form of a low-dimensional manifold in some highdimensional space dictated by the measurement apparatus (see Fig. 1). These so-called manifold-embedding approaches are in essence nonlinear (kernel) principal components techniques [26], which seek to display in some low-dimensional Euclidean space the manifold representing the correlations in the data. While powerful, such approaches face three challenges: computational cost; robustness against noise; and the assignment of physical meaning to the outcome of the analysis, i.e., physically correct interpretation.
Bayesian manifold approaches [19,27] and their equivalents [20,28] are able to operate at extremely low signal, but require prior knowledge of the manifold dimensional-ity and its physical meaning. They also display extremely unfavorable scaling behaviors [19][20][21]28]. It has thus not been possible, for example, to reconstruct objects with diameters exceeding eight times the spatial resolution, severely limiting the size of amenable objects and/or the resolution of the reconstruction.
Non-Bayesian graph-theoretic methods are computationally efficient, but tend not to be robust against noise [29,30]. More fundamentally, it has proved difficult to assign physical meaning to the outcome of graph-theoretic analyses [9], in the sense that the meaning of the dimensions in which the data manifold is embedded is unknown. Strategies to overcome this problem have included exploring the graph-theoretic consequences of specific changes in model systems [10], making restrictive assumptions about the nature of the data, and/or extracting specific information from the data first and subjecting this information to graph-theoretic analysis [9,11].
In this paper, we present a computationally efficient theoretical framework capable of interpreting the outcome of graph-theoretic analysis of scattering data without restrictive assumptions. Using this approach, we demonstrate 3D structure recovery from 2D diffraction snapshots of unknown orientation at computational complexities four orders of magnitude higher than hitherto possible [19,20]. In Paper II, we show that this framework can be used to: (1) recover structure from simulated and experimental snapshots at signal levels ∼ 10× lower than currently in use; and (2) reconstruct timeseries (movies) from ultralow-signal random sightings of evolving objects. In sum, therefore, our approach offers a powerful route to recovering structure and dynamics (3D movies) from ultralow-signal snapshots with no orientational or timing information.

IV. SYMMETRY-BASED ANALYSIS OF SCATTERING DATA
In this section we develop a mathematical framework for analyzing ensembles of 2D snapshots, using far-field scattering by a single object as a model problem to focus the discussion (see Sec. IV A). The basic principle of our approach, laid out in Sec. IV B, is that the data acquisition process can be described as a manifold embedding Φ [31,32], mapping the set of orientations of the object, SO(3), to the Hilbert space of snapshots on the detector plane. As a result, the differential-geometric properties of the rotation group formally carry over to the scattering dataset. In particular, the dataset can be equipped with a homogeneous Riemannian metric B, whose Laplacian eigenfunctions are the well-known Wigner D-functions [33][34][35]. Here, a metric is called homogeneous if any two points on the data manifold can be connected via a transformation that leaves the metric invariant. We refer to this class of transformations as symmetries. In Sec. IV C, we show that taking advantage of symmetry, as manifested in the properties of homogeneous metrics on SO(3) [36], leads to a powerful means for recovering snapshot orientations and hence the 3D structure of objects.
A number of sparse algorithms [5,8] are able to compute discrete approximations of Laplacian eigenfunctions directly from the data. However, these algorithms do not provide the eigenfunctions associated with B, but rather the eigenfunctions of an induced metric g associated with the embedding Φ (i.e., the measurement process). The properties of this induced metric are discussed in Sec. IV D. There, we show that g is not homogeneous, but admits a decomposition into a homogeneous metric plus a residual. Intriguingly, the homogeneous part of g corresponds to a well-known solution of general relativity (the so-called Taub solution [37]), which has the important property of preserving the Wigner D-functions as solutions of the Laplacian eigenvalue problem [36]. As described in Sec. IV E (and demonstrated in Sec. V and Paper II), this property applies to a broad range of scattering modalities, and can be exploited to perform highly accurate 3D reconstruction in a computationally-efficient manner.

A. Image formation
We first treat image formation as elastic, kinematic scattering of unpolarized radiation onto a far-field detector in reciprocal space [38], where each incident photon can scatter from the object at most once (kinematic scattering), and energy is conserved (elastic scattering). We show later that our conclusions are more generally applicable. In this minimal model, illustrated in Fig. 2, an incident beam of radiation with wavevector q 1 = Qz (we set Q > 0 by convention) is scattered by an object of density ρ(u) with Patterson function P (u) = du ρ(u )ρ(u − u ) = P (−u), where u and u are position vectors in R 3 . The structure-factor amplitude at point r on a detector plane fixed at right angles to the incident beam is given by the usual integral, a( r) = ω( r) du P (u) exp(iu · q( r)) where ω( r) is an obliquity factor proportional to the solid angle subtended at u = 0 by the detector element at r, and q( r) is the change in wavevector due to scattering. In elastic scattering we have where q 2 ( r) is the scattered wavevector, and θ and φ are the polar and azimuthal angles in reciprocal space corresponding to position r on the detector plane (see Fig. 2). Note that, by the convolution theorem, the Fourier transform F(P ) of the Patterson function is equal to the nonnegative function |F(ρ)| 2 . As a result, no modulus sign is needed in Eq. (1), and a 2 ( r j ) is equal to the intensity I j at pixel j in Fig. 1.
In an idealized, noise-free experiment involving a single object conformation, one observes a sequence of s snapshots on a detector of infinite extent, with the snapshots arising from random orientations of the object. Each snapshot is obtained from Eq. (1) by replacing P (u) with P R (u) = P (R −1 u), where R is a 3 × 3 right-handed rotation matrix; i.e., R T R = I and det(R) = 1.
The set of all matrices satisfying these conditions form the 3D rotation group, SO(3). It is well known that SO(3) is a Lie group, i.e., it is a differentiable manifold (in this case of dimension 3) [31,35,39]. Among the several parameterizations (coordinate charts) of SO(3), of interest to us here will be Euler angles, unit quaternions, and hyperspherical coordinates [40]. The last coordinate system stems from the fact that SO(3) has the topology of a three-sphere with its antipodal points identified.
For the remainder of this section, SO(3) will play the role of the latent manifold S, i.e., the set of degrees of freedom available to the object. In more general applications, the latent manifold would be augmented to contain the additional degrees of freedom, provided, of course, that these degrees of freedom admit a manifold description-a natural requirement for operations such as shifts and smooth conformational changes.

B. Riemannian formulation of image formation
We are concerned with intensity patterns generated by scattering, i.e., a subset of all possible patterns on a planar detector, referred to as data space. It is reasonable to expect that in a physical experiment involving a finitepower beam, the resulting distributions of structurefactor moduli belong to the set of square-integrable functions on the detector plane, L 2 (R 2 ). This is a Hilbert space of scalar functions f i equipped with the usual in-ner product and corresponding norm In many respects, L 2 (R 2 ) can be though of as a generalization of Euclidean space to infinite dimensions. In particular, the L 2 norm induces a distance f 1 − f 2 analogous to the standard distance in finite-dimensional Euclidean space. Moreover, viewed as a manifold [41], L 2 (R 2 ) has the important property that its elements are in one-to-one correspondence with its tangent spaces. Thus, the inner product in Eq. (3) can be interpreted as a metric tensor acting on pairs of tangent vectors on L 2 (R 2 ), or manifolds embedded in L 2 (R 2 ). We describe the image formation process as an embedding [31,32], Φ : S → L 2 (R 2 ), taking the latent manifold into data space. For instance, in the present application with S = SO(3), we have the explicit formula Φ(R) = a R , where R is an SO(3) rotation matrix and a R is the pattern on the detector given by Eq. (1) with P (u) replaced by P (R −1 u).
The image M = Φ(S) of the latent manifold in data space is called the data manifold (see Fig. 1). In the absence of degeneracies such as object symmetry, the data and latent manifolds are diffeomorphic manifolds, i.e., completely equivalent from the point of view of differential geometry. In contrast with manifolds of shapes [24,42] which can contain singularities, these are manifolds of operations and thus generally well-behaved. The image of the latent manifold in data space is then a smooth (here, three-dimensional) embedded hypersurface, which preserves the topology of the latent manifold, and the structure of its tangent spaces [43]. Perception, at least in this simple model, can be viewed as understanding the map between the latent manifold of orientations and the data manifold of pixel intensities.
The inner product induced on the data manifold by the inner product in data space leads to a metric tensor g on the latent manifold, which encodes the properties of the object and the imaging process. This metric is is constructed by converting the inner product of L 2 (R 2 ) in Eq. (3) to an equivalent inner product between tangent vectors on the data manifold. Specifically, given tangent vectors v 1 and v 2 on S, the induced metric g is defined through the action where Φ * is the so-called derivative map associated with Φ [41,44], carrying along tangent vectors on S to tangent vectors on M. The induced metric can be expanded in a suitable tensorial basis for S as a 3×3 symmetric positivedefinite (SPD) matrix with components g µν , viz., where {E 1 , E 2 , E 3 } is a basis of dual vector fields on S. Note that ds 2 = g(δv, δv) corresponds to the squared distance between two nearby points on S with relative separation δv. The integral of ds 2 over curves on S provides the data manifold with a notion of length and distance. For the purposes of the symmetry analysis ahead, it is useful to consider expansions of g in a basis of right-invariant vector fields [31,39], where, following the derivation in Appendix A, the components of g at orientation R are found to be In the above, dr denotes integration over the detector plane; q( r) is the change in wavevector given by Eq. (2); ∇ = (∂/∂q x , ∂/∂q y , ∂/∂q z ) is the gradient operator in reciprocal space; and J µ are the 3 × 3 antisymmetric matrices in Eq. (A5) generating rotations about the x, y, and z axes, respectively. As stated above, a R is the structure-factor amplitude corresponding to orientation R.

C. Symmetries
We now show that the Riemannian formulation described above reveals important symmetries, which can be used to determine the object orientation separately from the object itself. A fundamental concept in the discussion below is the notion of an isometry [31,44]. Broadly speaking, an isometry is a continuous invertible transformation that leaves g invariant. More specifically, any diffeomorphism φ : S → S mapping the latent manifold to itself induces a transformation φ * acting on tensors of the manifold; the latter will be an isometry if φ * (g) = g holds everywhere on S. A symmetry, therefore, in this context is an operation that leaves distances on the data manifold unchanged.
If the group of isometries of g acts transitively on S (i.e., any two points in S can be connected via a φ transformation), the pair (S, g) becomes a Riemannian homogeneous space [39]. We refer to any g meeting this condition as a homogeneous metric. Riemannian homogeneous spaces possess natural sets of orthonormal basis functions, analogous to the Fourier functions on the line and the spherical harmonics on the sphere, which can be employed for efficient data analysis [35]. It is therefore reasonable to design algorithms that explicitly take into account the underlying Riemannian symmetries of scattering data sets.
As discussed in more detail below, the isometry group of g in Eq. (5) is generally not large-enough to induce a transitive action; i.e., there exist points on the manifold that cannot be mapped to one another through an isometry. Nevertheless, considerable progress in the interpretation of datasets produced by scattering can be made by establishing the existence of a related metric h, whose isometry group meets the conditions for transitivity. As shown below, the properties of h can be used to interpret the results of data analysis performed on (S, g). By combining aspects of group theory and differential geometry, the general techniques developed here constitute a novel approach for analyzing scattering data, and potentially other machine-learning applications.
The canonical classes of symmetry operations in problems involving orientational degrees of freedom are the so-called left and right multiplication maps [31,39,44]. These maps, respectively denoted L Q and R Q , are parameterized by an arbitrary rotation matrix Q in SO (3), and act on SO(3) elements by multiplication from the left or the right, respectively. That is, Viewed as groups, the collection of all left multiplication maps or right multiplication maps are isomorphic to SO (3), from which it immediately follows that their action is transitive. Thus, any metric g invariant under L Q or R Q (or both) makes (S, g) a Riemannian homogeneous space.
The natural metric for the SO(3) latent manifold of orientations with these symmetries is the metric tensor B associated with the Killing form of the Lie algebra of SO(3) [31,39]. This metric is bi-invariant under left and right translations. The corresponding group of isometries has SO(3) × SO(3) structure, i.e., it is a six-dimensional group. In hyperspherical coordinates (χ, θ, φ), B is represented by the diagonal matrix which is identical to the canonical metric on the threesphere S 3 . For this reason B is frequently referred to as the "round" metric on SO (3).
A key implication of the symmetries of B pertains to the eigenfunctions of the corresponding Laplace-Beltrami operator ∆ B [45], defined as with positive integer j and integers m and m in the range [−j, j]. Written out in terms of Euler angles in the zyz convention, (α 1 , α 2 , α 3 ), explicit formulas for the j = 1 D-functions are As one can check by algebraic manipulation, the above ninefold-degenerate set of eigenfunctions can be put into one-to-one correspondence with the elements of R. That is, it is possible to find linear combination coefficients c pqmm such that where R pq are the elements of R. Thus, if eigenfunctions of ∆ B could be accessed by processing the observed snapshots a R on the data manifold, e.g., through one of the graph-theoretic algorithms developed in machine learning [5,8], Eq. (13) could be used to invert the embedding map Φ. This would be tantamount to having learned to "navigate" on the data manifold. The method outlined above for symmetry-based inversion of Φ is also applicable under weaker symmetry conditions on the metric. In particular, it can be extended to certain metrics of the form where E µ are the right-invariant dual basis vectors in Eq. (A8), and µ are non-negative parameters. These are the so-called spinning-top metrics of classical and quantum-mechanical rotors [34,48], arising also in homogeneous cosmological models of general relativity [36,37,49]. In classical and quantum mechanics, h features in the Hamiltonian of a rotating rigid body with principal moments of inertia given by µ . In general relativity, µ characterize the anisotropy of space in the so-called mixmaster cosmological model [49]. By construction, metrics in this family are invariant under arbitrary right multiplications, i.e., they possess an SO(3) isometry group. This is sufficient to make (S, h) a Riemannian homogeneous space. If two of the µ parameters are equal (e.g., 1 = 2 ), h describes the motion of an axisymmetric rotor, or the spatial structure of the Taub solution in general relativity [37]. As one may explicitly verify, in addition to being invariant under left multiplications, the metric of the axisymmetric rotor is invariant under translations from the left by the oneparameter subgroup associated with rotations about the axis of symmetry, which corresponds here to the direction of the incoming scattering beam. Therefore, it has a larger, SO(3) × SO(2) isometry group than the more general metric of an asymmetric rotor. In the special case that all of the µ are equal, h reduces to a multiple of the round metric B in Eq. (9).
Let ∆ h denote the Laplace-Beltrami operator associated with h. It is possible to verify that ∆ h , and the Laplacian ∆ B corresponding to the round metric commute. That is, it is possible to find eigenfunctions that satisfy the eigenvalue problem for ∆ B and ∆ h simultaneously. In particular, as Hu [36] shows, the eigenvalueeigenfunction pairs (λ j m , y j m ) of ∆ h can be expressed as linear combinations of Wigner D-functions with the same j and m quantum numbers: for some linear expansion coefficients A j mm . For the most general metrics with 1 = 2 = 3 , no closed-form expression is available for either A j mm , or the corresponding eigenvalues λ j m . However, in the special case of axisymmetric rotors, any A j mm in Eq. (15) will produce an eigenfunction of ∆ h . This means that the vector space spanned by the j = 1 eigenfunctions associated with the metric of an axisymmetric rotor (denoted here by y j mm ) is nine-dimensional. In particular, Eq. (13) can be applied directly with D 1 mm replaced by y 1 mm . The eigenvalue λ j m corresponding to y j mm in the socalled prolate configuration, 1 = 2 ≤ 3 , is given by with a similar formula for the oblate configuration 1 = 2 ≥ 3 [36]. An important point about Eq. (16) (and its oblate analog) is that all of the j = 1 eigenvalues are greater than zero and smaller than the maximal j = 2 eigenvalue; i.e., the j = 1 eigenfunctions can be identified by ordering the eigenvalues. This property, which does not apply for general asymmetric-rotor metrics, alleviates the difficulty of identifying the correct subset of eigenfunctions for orientation recovery via Eq. (13).

D. Accessing the homogeneous metric
The method outlined above for symmetry-based inversion of Φ makes direct use of the fact that (S, h) is a Riemannian homogeneous space. However, (S, h) is generally inaccessible to numerical algorithms, which access a discrete subset of the latent manifold equipped with the induced metric g, i.e., (S, g). It is therefore important to establish the isometries of g, because, for instance, the latter govern the behavior of the Laplacian eigenfunctions computed via graph-theoretic algorithms [3-5, 7, 8].
In general, g is not invariant under arbitrary left or right multiplications. The lack of invariance of g under the right multiplication map R Q can be seen by considering its components in a right-invariant basis. In particular, it can be shown that R Q is an isometry of g if and only if the components g µν in the right-invariant basis from Eq. (7) are invariant under replacing R by RQ, which does not hold generally.
What about the behavior of g under left multiplications? Here, an expansion of g in an analogous leftinvariant basis reveals that g is not invariant under left multiplications by general SO(3) matrices, but is invariant under left multiplication by a rotation matrix Q about the beam axis z (see Fig. 2). This is intuitively obvious, because the outcome of replacing a snapshot a R with a QR is then a rotation on the detector plane, which has no influence on the value of the L 2 inner product used to evaluate g. In summary, instead of the six-dimensional SO(3) × SO(3) isometry group of B, the isometry group of the induced metric is the one-dimensional SO(2) group of rotations about the beam axis, which obviously cannot act transitively on the three-dimensional latent manifold.
Despite the insufficiency of g to make (S, g) a Riemannian homogeneous space, g admits a decomposition of the form with the following key properties: (1) h is a spinning-top metric from Eq. (14) with 1 = 2 (i.e., an axisymmetricrotor metric, not an arbitrary right-invariant metric). (2) w is a symmetric (but not necessarily positive-definite) tensor, whose components average to zero over the manifold, describing the inhomogeneous part of g (see Appendix A 2 for a derivation). Moreover, provided that w meets a suitable norm constraint, g and h are bi-Lipschitz equivalent. This means that the distortion in distances on S caused by using g instead of h, and the corresponding change in the Laplace-Beltrami eigenfunctions, are both bounded [50]. The fact that h in Eq. (17) is the metric of an axisymmetric rotor is a direct consequence of projection onto a circularly-symmetric 2D detector, and therefore an essential aspect of scattering experiments. The important point is this: scattering manifolds are associated with induced metrics g, which, in themselves, possess low symmetry. But they can be decomposed into a homogeneous part h with a high SO(2) × SO(3) symmetry and thus associated with well-known Laplacian eigenfunctions independently of the object, plus a low-symmetry residual, which depends on the object. In essence, one can regard the homogeneous metric h as an object-independent "Platonic ideal" [2] version of g, and develop algorithms, which make use of the properties of h to analyze all scattering datasets.
The symmetries of h can be exploited in several ways: (1) enforce appropriate constraints in Bayesian algorithms [19]; (2) interpret the results of numerical Ricci flow [51,52]; or (3) approximate the eigenfunctions associated with the homogeneous metric by a truncated set of eigenfunctions associated with the inhomogeneous metric. The last possibility is described below.
Represent each snapshot a Ri in the dataset as a point in nine-dimensional Euclidean space R 9 with coordinates given by (ψ i1 , . . . , ψ i9 ). Here, ψ ik are eigenvector components of an elliptic operator on a graph constructed from the s observed snapshots (a R1 , . . . , a Rs ) in n-dimensional data space (see Appendix B). This induces an embedding Ψ of the latent manifold S in R 9 given by Provided that the number of observed snapshots is sufficiently large, suitable algorithms (such as Diffusion Map of Coifman and Lafon [8,53]), lead to embedding coordinates ψ ik , which converge to the corresponding values of the Laplacian eigenfunctions associated with g, even if the sampling of the data manifold is non-uniform. More specifically, for large-enough s we have the correspondence where ψ k (R i ) is the k-th eigenfunction of the Laplace-Beltrami operator ∆ g associated with the induced metric in Eq. (5); i.e., Now we explicitly employ the symmetries of the homogeneous metric h in Eq. (14), motivated by the decomposition of g in Eq. (17): Express each eigenfunction of the homogeneous Laplacian ∆ h with j = 1 as a linear combination of the eigenfunctions ψ k of the inhomogeneous Laplacian, viz., where we have made use of the correspondence in Eq.

E. Extension to general scattering problems
Because the above analysis was performed under a restrictive set of experimental conditions, we now consider the effect of removing these. The assumption of kinematic (single) scattering introduces an additional symmetry due to Friedel's law. As we have not used this symmetry, our arguments remain valid under multiple scattering conditions. The use of linearly or elliptically polarized radiation introduces a second preferred direction in addition to that of the beam, thus removing the SO(2) isometry under left translations, but this can be restored by appropriate correction with the polarization factor. A detector not at right angles to the beam axis also reduces the isometry to I × SO(3), but this can also be easily corrected by an appropriate geometric factor. Absorption can be accommodated as a complex density function, inelastic scattering by allowing q 1 = q 2 , etc. Neither affects our conclusions. Our approach is thus applicable to a wide range of image formation modalities.
Note, however, that the method has to be modified when the object has discrete or continuous symmetries. This is because the data manifold in this case is not SO(3), but the quotient space SO(3)/Γ where Γ is a subgroup of SO(3) representing the object's symmetries. Among the eigenfunctions of the Laplacian on the SO(3) manifold, only those that are constant on Γ "survive" in the SO(3)/Γ environment. Thus, the orientation recovery procedure must be modified depending on the form of the available eigenfunctions. As mentioned in Sec. II, this issue will be addressed elsewhere.

V. DEMONSTRATION OF STRUCTURE RECOVERY IN 3D
It has long been known that the use of problem-specific constraints can substantially increase computational efficiency [54]. By combining wide applicability with class specificity, symmetries represent a particularly powerful example of such constraints. Exploiting these, we here demonstrate successful orientation recovery for a system computationally 10 4 × more complex than previously attempted [19,20]. In Paper II we apply our framework to snapshots of various kinds with extremely low signal.

A. 3D reconstruction from diffraction snapshots with high computational complexity
We simulated X-ray snapshots of the closed conformation of E. coli adenylate kinase (ADK; PDB descriptor 1ANK) in different orientations to a spatial resolution d = 2.45Å, using 1Å photons. In this calculation, Cromer-Mann atomic scattering factors were used for the 1656 non-hydrogen atoms [55], neglecting the hydrogen atoms. We discretized the diffraction patterns on a uniform grid of n = 126 × 126 = 15, 876 detector pixels of appropriate (Shannon) size, taking the corresponding orientations uniformly on SO(3) according to the algorithm in Ref. [56]. The number of diffraction patterns required to sample SO(3) adequately is given by 8π 2 (D/d) 3 ≈ 8.5 × 10 5 , where D = 54Å is the diameter of the molecule [19]. In our simulation, however, a larger D = 72Å diameter was assumed, allowing, e.g., for the possibility to reconstruct the structure of ADK's open conformation. Hence, a total of s = 2×10 6 patterns were used in the present analysis.
The diffraction patterns were provided to our Diffusion Map-based algorithm with no orientational information, and the orientation of each diffraction pattern was determined by means of the algorithm in Table III [57]. To estimate the difference between the deduced and true orientations, respectively represented by unit quaternionsτ i and τ i [see Eq. (B9)], we used the RMS internal angular distance error, In the above, D ij = 2 arccos(|τ i · τ j |) andD ij = 2 arccos(|τ i ·τ j |) are respectively the true and estimated internal distances between orientations R i and R j , and · is the inner product between quaternions. The resulting internal angular distance error in the present calculation was 0.8 Shannon angles. Next, we placed the diffracted intensities onto a uniform 3D Cartesian grid by "cone-gridding" [58], and deduced the 3D electron density by iterative phasing with a variant of the charge flipping technique [16] developed by Marchesini [59]. Charge flipping was prevented outside a spherical support about twice the diameter of the molecule. The R-factor between the gridded scattering amplitudes and the ones obtained from phasing was 0.19. The close agreement with the correct structural model is shown in Fig. 3 and the EPAPS movie.

B. Computational cost
In the absence of orientational information, the computational cost C of orientation recovery without restrictive sparsity assumptions scales as a power law with D the object diameter and d the spatial resolution [19][20][21]. The analysis of ADK was performed with R = 30, compared with the largest previously published value of R ≤ 8 [19][20][21]. This represents an increase of four orders of magnitude in computational complexity over the state of the art, as shown below. The computational cost of a single expectationmaximization (EM) step in Bayesian algorithms (e.g., the GTM algorithm [27,60]) scales as Ksn, where K, s, and n are the number of quaternion nodes, the number of snapshots, and the data-space dimension, respectively [19]. Here, the data space dimension is equal to the number of Shannon pixels. Assuming an oversampling of 2 in each linear dimension, we have where S is the object symmetry, f the number of snapshots per orientational bin, and λ the wavelength. Assuming that the number of EM iterations is constant, Fung et al. [19] estimate an eighth-power scaling of the form in Eq. (24) for the total computational cost of orientation recovery with GTM. The fact that GTM is NPhard remains moot. Loh and Elser [20] argue that their EM-based approach scales as where M rot is the number of rotation samples (orientational bins), and that this leads to an R 6 scaling. This is based on three assumptions: (1) A sparse formulation can be used to replace the number of detector pixels n with a much smaller number of scattered photons N photons . (2) N photons does not depend on R, and thus can be ignored in the scaling behavior. (3) The number of EM iterations is independent of R (as in Ref. [19]), and of N photons . Note that M rot and the number of GTM bins, K, are equal. Assumption (1) holds only for very small photon counts and in the absence of significant background scattering. Assumption (2) is not justified [61]. For a globular protein, the total number of photons scattered to large angle scales linearly with the number of non-H atoms, and hence object volume, i.e., as D 3 . Thus, (28) Using the scaling expression of Eq. (26), the increase in complexity on going from GroEL to ADK (72Å support) for the Loh and Elser algorithm varies between O(10 3 ) and O(10 4 ) for sparse and non-sparse representations, respectively (see Table I). As a sparse representation is not generally possible and we have not had to resort to it, the appropriate comparison is 10 4 .

VI. CONCLUSIONS
We have shown that a Riemannian formulation of scattering reveals underlying, object-independent sym-metries, which stem from the nature of operations in 3D space and projection onto a 2D detector. These symmetries lead to immediate identification of the natural eigenfunctions of manifolds produced by data from a wide range of scattering scenarios, and thus to physicallybased interpretation of the outcome of graph-theoretic analysis of such data. In practical terms, the ability to access the homogeneous metric offers a computationally efficient route to determining object orientation without object recognition, while the object-dependent term provides a concise fingerprint of the object for recognition purposes. There are tantalizing indications that face perception in higher primates may occur in this way [1]. The ability to use symmetries to navigate on the data manifold is tantamount to efficient machine learning in 3D, in the sense that given any 2D projection, any other can be reconstructed.
As shown in Paper II, the manifold itself offers a powerful route to image reconstruction at extremely low signal, because snapshots reconstructed from the manifold achieve higher signal-to-noise ratios than can be obtained by traditional methods relying on classification and averaging. Combined with the ability to sort random snapshots of an evolving system into a time-series, also demonstrated in Paper II, our approach offers a radically new way for studying dynamic systems in 4D. Fundamentally, the homogeneous metric describes the transformations of objects without reference to any specific object. This is reminiscent of a Platonic Form, from which specific objects emanate [2]. It is therefore tempting to regard the homogeneous manifold as a Platonic Form, from which our perception of three-dimensional objects stems [62]. Here, we discuss the properties of the induced metric tensor g in Eq. (5). We begin in Appendix A 1 with a derivation of the explicit form for the components of g in the right-invariant basis from Eq. (7). In Appendix A 2 we perform the decomposition of g in Eq. (17) into axisymmetric and inhomogeneous parts. In Appendix A 3 we provide estimates of the error incurred in using the Laplacian eigenfunctions associated with the inhomogeneous metric g for orientation recovery.

Derivation
Following Sec. IV B, we describe scattering as an embedding Φ taking the latent manifold S to the set of square-integrable intensity patterns L 2 (R 2 ) on a 2D dimensional detector. Broadly speaking, the induced metric describes an inner product between tangent vectors on the latent manifold S associated with that embedding. More specifically, in accordance with Eq. (5), that inner product is computed by "pushing forward" tangent vectors on S to manifest space, and applying the canonical Hilbert space inner product in Eq. (3). The map carrying along the tangent vectors of S is the derivative map Φ * associated with Φ, which is evaluated as follows.
First, note that every smooth tangent vector field v on S generates a corresponding one-parameter family of transformations [31,32,39], with α a scalar parameter and R α an element of S, which depends smoothly on α. In the above, φ α describes a curve on the manifold (called an integral curve of v), which is tangent to v at every point. Likewise, the pushforward Φ * (v) of v is associated with a continuous transformationφ α of intensity snapshots in manifest space. Denoting the snapshot associated with orientation R by a R = Φ(R), that transformation is given byφ with a Rα determined from Eq. (A1). The outcome of acting on v with Φ * is then the directional derivative of intensity snapshots along the path defined byφ α , viz.
The induced metric g resulting from the above procedure will depend on the explicit form of the embedding map Φ. Here, we consider the case where Φ describes far-field kinematic elastic scattering from a single object, where intensity amplitudes are given by the Fourier integral in Eq. (1). Other scattering scenarios can be treated in a similar manner, provided that Φ meets the conditions of an embedding. For our purposes, it is sufficient to consider one-parameter transformations arising from left multiplication by SO(3) matrices carrying out rotations about one of the x, y, or z axes [see Eq. (8)]. That is, we set where α is the rotation angle in radians, and J 1 , J 2 , and J 3 are 3 × 3 antisymmetric matrices generating rotations about the x, y, and z axes: We denote the vector fields generating φ µ α by e µ . It then follows from Eq. (A3) with R α = L µ α (R) that the pushforward fields Φ * (e µ ) are given by Note that in deriving Eq. (A6) we have used the divergence-free property ∇ · (J µ q) = 0, which applies for any scattering wavevector q and rotation generator J µ . As one may verify, the tangent vector fields e µ are linearly independent. This means that the set {e 1 , e 2 , e 3 } forms a basis to expand tangent vectors on S, and, in turn, the induced metric can be represented by a 3 × 3 matrix with elements g µν (R) = (Φ * e µ , Φ * e ν ). (A7) The g µν (R) above are the components of the induced metric in Eq. (6) at orientation R, provided that E µ are dual basis vectors to e ν , defined through the relation Moreover, it is possible to show that e µ are invariant under the right multiplication map R Q in Eq. (8) [44]. Substituting for Φ * (e µ ) in Eq. (A7) using Eq. (A6) then leads to the expression in Eq. (7) for the components of the induced metric in a right-invariant tensorial basis. Besides the form in Eq. (7), it is useful to express the components of g in terms of spherical harmonics in reciprocal space. As is customary in diffraction theory [63], we write down the reference scattering amplitude a = Φ(I) corresponding to the identity matrix I using the expansion where a m j are complex functions of the radial coordinate q in reciprocal space, and Y m j are spherical harmonics. Note the property a m * j = (−1) m a −m j , which is a consequence of a(q) being real. Making use of the standard formula describing rotations of spherical harmonics via Wigner D-matrices, it follows that the amplitude distribution a R = Φ(R) corresponding to object orientation R is given by Inserting the above in Eq. (7) leads to the following general expression for the metric components, In the above, K

Decomposition into homogeneous and inhomogeneous parts
Even though the induced metric tensor g in Eq. (5) is not homogeneous, it is nevertheless possible to decompose it as a sum where h is a homogeneous metric with respect to a group acting transitively on the data manifold, and w a tensor that averages to zero over the manifold. Using Fourier theory on SO(3) [35,64], here we compute the components of h and w in the right-invariant basis E µ from Eq. (A8). The analysis presented below yields the important results that (i) the homogeneous part h of the metric belongs to the family of metrics associated with axisymmetric rotors; (ii) the Fourier spectrum of w has a highly sparse structure, with only a limited number of nonzero coefficients contributing to metric inhomogeneity. If, as a result of the structure of its Fourier spectrum, the norm w of w (suitably defined) is smaller than the norm h of h, then using the numerically-accessible eigenfunctions of g to approximate the eigenfunctions of h results in a bounded loss of accuracy.
We begin by noting that the Fourier transform of the metric components g µν in Eq. (A12) consists for each (µ, ν) of the sequence {ĝ j µν } ∞ j=0 of (2j + 1) × (2j + 1) matrices g j µν , whose components are given by In the above, integration is performed over SO(3) using the volume element dV of the bi-invariant metric B in Eq. (9) [65], and D j mm are the Wigner D-functions from Eq. (11). The collection of theĝ j µν matrices may be used to recover g (R) via the inverse transform where D j (R) = [D h mm (R)] is the j-th Wigner matrix of size (2j + 1) × (2j + 1). Because D 0 is the trivial representation of SO(3), a constant scalar on the manifold, the termĝ 0 µν corresponds to the homogeneous part of the metric. Theĝ j µν matrices with j ≥ 1 give rise to the inhomogeneous tensor w, which encodes object-specific information. That is, we have (2j + 1) tr(ĝ j µν D j (R)).
(A17) where h µν and w µν are the components of h and w in the E µ basis, respectively.
Since each component of the metric contains a product of two Wigner D-functions [see Eq. (A12)], the calculation of the Fourier transform is facilitated significantly by the triple-integral formula involving the 3-jm symbols of quantum mechanical angular momentum [35,47], In particular, inserting the metric components from Eq. (A12) in Eq. (A15), and making use of the triple-integral formula, leads to the result where we have employed the triangle inequality, |j 1 − j 2 | ≤ j ≤ j 1 + j 2 , obeyed by the nonzero 3-jm symbols.
The (A20) It follows that the j = 0 term of the Fourier spectrum (defined only for m = m = 0) giving rise to the homogeneous part of the metric has no off-diagonal components. This in turn means that h can be expressed in terms of some non-negative parameters µ in the form of Eq. (14). An explicit evaluation of the on-diagonal terms yields the additional result that the 1 and 2 parameters are, in fact, equal. Thus, the components of the induced metric read By construction, the average of w µν over the manifold vanishes, i.e., dV (R) w µν (R) = 0. (A22)

Bounding the error in the eigenfunctions
Given the above decomposition of the induced metric g into a homogeneous metric h associated with an axisymmetric rotor plus an inhomogeneous component w, it is possible to write down sufficient conditions for the discrepancy between the Laplace-Beltrami eigenfunctions associated with g and h to be bounded. First, note that h induces a norm · h,∞ on tensors of type (0, 2) on the latent manifold. For a sufficiently smooth tensor field u of type (0, 2) that norm is defined as where the supremum is taken over SO(3) matrices R and nonzero tangent vectors v on SO(3). Thus, u h,∞ involves computing the least upper bound over the manifold and over the nonzero tangent vectors of the relative magnitudes of the quadratic forms u(v, v) and h(v, v). Next, set u = g = h + w, and compute By the direct and reverse triangle inequalities of norms, we have the bounds If w h,∞ is strictly less than one, then Eq. (A25) implies that g h,∞ is bounded from below by the nonzero constant Λ − = 1 − w h,∞ and from above by It follows that if the condition w h,∞ < 1 holds, then g and h obey the bound This property-is called bi-Lipschitz equivalence. Bérard et al. [50] show that if the Ricci curvatures of two bi-Lipschitz equivalent metrics, g and h, are bounded from below by some negative constant, then there exist constants η i (Λ), which go to zero as Λ approaches 1, such that for any orthonormal basis of eigenfunctions {y i } K i=1 of ∆ h , one can associate an orthonormal basis of eigenfunctions {ψ i } K i=1 of ∆ g satisfying the bound where · ∞ denotes the supremum norm. This means that if one approximates the set of (possibly degenerate) eigenfunctions of ∆ h with corresponding eigenvalues λ h 1 , . . . , λ h K by the eigenfunctions of ∆ g with corresponding eigenvalues λ 1 , . . . , λ K (sorted in increasing order), then the maximum error is bounded from above in a pointwise sense. Here, we do not attempt to derive conditions on the properties of the object needed to guarantee that w h,∞ < 1 holds, nor do we address whether the bi-Lipschitz equivalence between g and h can be deduced by imposing weaker constraints on w. Nevertheless, to the extent that the norm of w is indeed small (as suggested by, e.g., the strong selection rules on its Fourier spectrum), we expect that using the leading eigenfunctions of g for orientation recovery should produce only small systematic error, as observed in practice.

Appendix B: Algorithms
As an instructive application of the symmetries identified in Sec. IV, we describe an accurate and efficient algorithm for the analysis of scattering datasets. This makes explicit use of the properties of the homogeneous metric h in Eq. (14) to interpret results produced by the Diffusion Map algorithm of Coifman and Lafon [8]. As mentioned earlier, and also demonstrated in Paper II, a variety of other manifold algorithms can also be used, all taking advantage of the symmetries underlying datasets produced by scattering. Here, we consider the idealized scenario of noise-free data. This will be relaxed in Paper II, where we extend our approach to deal with snapshots severely affected by Poisson and Gaussian noise.
Instead of a continuous data manifold M, experimental data represent a countable subset M of M consisting of s identically and independently distributed (IID) samples in n-dimensional data space (with n the number of detector pixels) drawn randomly from a possibly non-uniform distribution on M. That is, we have where a i = (a i1 , . . . , a in ) are n-dimensional vectors of pixel amplitudes (see Fig. 1). As described in Sec. IV B, the amplitudes are given by a ij = a Ri ( r j ), where r j is the position of pixel j in the detector plane, and a Ri = Φ(R i ) is the snapshot associated with orientation R i . In Ref. [8] it is shown that it is possible to construct a one-parameter family of diffusion processes (random walks) on the point cloud from Eq. (B1), with each process described by an s × s transition probability matrix P , such that in the limit → 0 and s → ∞ the eigenvectors of P converge to the eigenfunctions of the Laplace-Beltrami operator ∆ g associated with the distance metric in data space [i.e., the induced metric tensor in Eq. (5)]. More specifically, if ψ k are s-dimensional column vectors in the eigenvalue problem then the relation ψ ik ≈ ψ k (R i ) holds for large-enough s and small-enough , where ψ k (R i ) is the k-th eigenfunction of ∆ g in Eq. (21) evaluated at element R i on the latent manifold.
Central to the efficiency of Diffusion Map is the fact that the transition probability matrix P is highly sparse. This is because P is constructed by suitable normalizations of a matrix W assigning weights W ij = K(a i , a j ) via a Gaussian kernel to pairs of snapshots in M , depending on their Euclidean distance in data space. In applications, one typically fixes a positive integer d and for each snapshot a i retains the distances up to its d-th nearest neighbor. Hereafter we denote the index of the j-th nearest neighbor of a i by N ij , and the corresponding distance by S ij = a i − a Nij . Pseudocode for computing P given the s × d distance and index matrices, S = [S ij ] and N = [N ij ], is listed in Table II.  Once the eigenvectors in Eq. (B2) have been evaluated, the next step in the orientation-recovery process is to evaluate the linear-combination coefficients needed to convert the nine-dimensional embedding coordinates (ψ l1 l, . . . , ψ l9 ) in Eq. (18) to elements of an approximate rotation matrixR l . This conversion is performed by the discrete analog of Eq. (13), namely [R l ] ij = 9 k=1 c ijk ψ lk l, 1 ≤ l ≤ s.
Here the expansion coefficients c ijk are determined by a least-squares minimization of the error functional where M = Because the eigenfunctions ψ k are not exact linear combinations of the Wigner D-functions in Eq. (12), a (real) polar decomposition step may be applied to produce an exact rotation matrix R, given an approximate rotation matrix estimated from the eigenfunction data via Eq. (B4). The real polar decomposition is a factorization of a real square matrixR into the product where R is an orthogonal matrix, and S is symmetric positive-semidefinite matrix.
The nonlinear least-squares step may be adequately performed using only a small subset of the full dataset, consisting of r datapoints drawn randomly from M . For instance, in the results of Sec. V, where the total number of samples is s = 2 × 10 6 , we use r = 8 × 10 4 datapoints. We denote the total squared residual of the optimization process by G * , where G * is equal to G({c ijk }) in Eq. (B5a), with {c ijk } given by the minimizer of the error functional.
Besides R, one might additionally require an evaluation of the corresponding coordinates in an SO(3) coordinate chart (e.g., for diffraction-pattern gridding). Depending on the coordinate chart of interest, various techniques are available to perform this procedure stably and efficiently. For instance, the coordinates of R in the unit-quaternion parameterization [40] can be determined by making use of the trace formula tr(R) = 1 + cos χ (B7) and the eigenvector relation where χ is the angle of the rotation represented by R, and n = n x x+n y y +n z z is a unit vector in R 3 directed along the corresponding axis of rotation. The unit quaternion τ corresponding to R is given by τ = (cos(χ/2), sin(χ/2)n x , sin(χ/2)n y , sin(χ/2)n z ). (B9) In outline, our orientation-recovery method consists of the following three basic steps: 1. Use Diffusion Map to compute the first nine nontrivial eigenfunctions ψ 1 , . . . , ψ 9 of the diffusion operator P ; 2. Determine the linear combination coefficients c ijk to transform the ψ k eigenvectors to (approximate) rotation matrices by solving the nonlinear least squares problem in Eq. (B5); 3. Use the correspondence in Eq. (B4) and polar decomposition [Eq. (B6)] to assign a rotation matrix R l to each observed diffraction pattern, and (optionally) extract the parameters of R l in the SO(3) parameterization of interest.
A high-level description of the orientation-recovery procedure for noise-free data is given in Table III. Before closing this section, we comment briefly on methods for choosing the kernel width and the number of retained nearest neighbors d. This is a common issue in manifold-embedding algorithms, and a number of strategies for determining these parameters have been developed in the literature [10,30]. Here, we determine and d a posteriori, by monitoring the value of the leastsquares residual G * . Specifically, in the applications presented in Sec. V, we start with a number of nearest neighbors d ∼ 100 and a value of of order the mean distance 2: return N 3: Compute the sparse transition probability matrix P using the algorithm in Table II with inputs S, N, , and α = 1.

11:
return τi 12: end for to the (d/2)-th nearest neighbors of the datapoints, and then perform successive refinements of these parameters seeking to minimize G * .