Three dimensional structure from intensity correlations

We develop the analysis of x-ray intensity correlations from dilute ensembles of identical particles in a number of ways. First, we show that the 3D particle structure can be determined if the particles can be aligned with respect to a single axis having a known angle with respect to the incident beam. Second, we clarify the phase problem in this setting and introduce a data reduction scheme that assesses the integrity of the data even before the particle reconstruction is attempted. Finally, we describe an algorithm that reconstructs intensity and particle density simultaneously, thereby making maximal use of the available constraints.


Introduction
As first pointed out by Kam [1], intensity fluctuations in an x-ray experiment with ensembles of non-oriented identical particles can provide structural information about the particles themselves. Correlations in the fluctuations arise from sets of photons that are scattered from the same particle and as such effectively provide a signal from individual particles. This technique does not place demands on the spatial coherence of the source -provided it spans at least one particle -but requires short pulses of high fluence and was largely not developed for that reason. With the availability of high intensity free-electron laser sources in recent years the situation has changed, and there is renewed interest in developing Kam's idea into a practical method. This effort has been led by Dilano and co-workers, with theoretical contributions [2,3] as well as proof-of-principle experiments [4,5].
There is currently no known method for extracting 3D structure from intensity fluctuation data. Moreover, a constraint counting argument [6] shows that information derived from experiments with completely non-oriented particles is deficient for extracting 3D structure. In this paper we revive the prospects of Kam's idea in the 3D realm provided the particles in the ensemble can be partially oriented. The type of alignment required is illustrated with pyramidal particles in Figure 1. As a result of native or engineered interactions and geometry, we propose that the particles adsorb on a substrate in such a way that a particular oriented axis within the particle is aligned with the substrate normal with near 100% probability 1 . Other means of achiev-Figure 1: Partially oriented particle ensemble realized as dice with a strong tendency to land on one particular face. Alternatively, partial alignment might be realized in the manner of membrane proteins embedded in a lipid membrane.
ing this degree of partial orientation would achieve the same goal. The concentration of their energy density makes substrates more efficient alignment mechanisms than laser fields, but this comes at the expense of introducing background and making certain ranges of the orientation angle inaccessible to the experiment. Membrane protein complexes, axially oriented in a synthetic lipid membrane, a liquid substrate, would be natural targets for this technique [3].
The fluctuations that form the basis of Kam's idea are not due to photon number fluctuations but arise from different realizations of the particle orientations. In our case, these are limited to rotations about the substrate normal. If the particles undergo rotational diffusion, then the duration of each intensity measurement must be short compared to the diffusion time scale. When diffusion is absent, the sample substrate must be translated in order that a different set of particles, with different orientations, intercepts the beam. In either case, and from a signal perspective, it is always advantageous to focus the beam on the smallest number of particles N , if this can be done so without adversely affecting beam divergence and the integrated fluence. That is because, while the mean intensities recorded at two positions on the detector are independent of N , their covariance -the source of structure information -is diminished by 1/N . Even when the background signal at the two pixels is perfectly uncorrelated, in being independent of N it will dominate the error in the covariance estimation for large N . There are other factors, however, that rule out small N even when this does not sacrifice the quality of the beam. Chief among these is particle damage, which must be kept at a low rate if subsequent measurements are taken on a single set of rotationally diffusing particles.
In a nutshell, our reconstruction method involves three stages. In the data collection stage, intensity correlations are recorded for five continuous parameters. Four of these are the x and y positions of two pixels on the detector and the fifth is the tilt angle θ of the substrate. This 5D data set is reduced to a 3D data set in the data reduction stage. The 3D data comprise complex numbers I m (r, z), or angular harmonics, giving the Fourier series expansion coefficients of the particle intensity on families of circles in a frame fixed to the particle. The data reduction also provides an internal consistency criterion, similar to principal component analysis, for establishing ranges for the radial (r) and longitudinal (z) frequency content, as well as the maximum harmonic order M that can be extracted from the data. The data reduction stage leaves M/2 phases in the intensity reconstruction undetermined. Finally, in the particle reconstruction stage, two phase problems are solved jointly that determine the 3D particle density from the (3D) angular harmonics. In addition to the usual phases encountered when reconstructing from intensity data, the M/2 unknown phases in the angular harmonics from the second stage are determined from the property that the intensity is consistent with a particle of compact support.

Data collection
The coordinate frame for data collection is defined by the incident x-ray beam, as z-axis, and orthogonal x and y axes in the plane of the detector. With the sample located at the origin, the unscattered beam would hit the detector at x = y = 0, z = Z. The axes on the detector are chosen with reference to the axis about which the substrate holding the sample is rotated. For tilt angle θ = 0 the substrate and detector planes are parallel and in-plane axes x ′ and y ′ on the substrate coincide with those on the detector. In all the measurements the substrate is tilted about the y ′ axis, leading to the following relationship between the lab and substrate axes: The z ′ axis is the substrate normal and coincides with an alignment axis fixed in the body of each particle adsorbed on the substrate. A third set of axes, fixed in a particular particle, is free to rotate with respect to the substrate about its z ′ axis by some angle α: x ′ = cos α x ′′ + sin α y ′′ (2a) In the experiment there is no control over the angle α and we assume it takes uniformly distributed random values for each particle.
We can now express the photon wave vector change, seen in the laboratory, in the frame fixed to a particular particle. For intensity correlation measurements we are interested in the intensities at two points on the detector plane, that we label 1 and 2. The source of the scattered photons is the same particle, since the covariance vanishes for different particles whose rotation angles α are assumed to be independent. Using equations (1) and (2), the wave vector transferred to photons detected at coordinates x 1 and y 1 on the detector is q 1 = (x 1 cos θ cos α − y 1 sin α) x ′′ + (x 1 cos θ sin α + y 1 cos α) y ′′ + (x 1 sin θ) z ′′ where λ is the wavelength and there is a similar equation for q 2 . In subsequent equations we will express all wave vectors, as here, in units of 2π/(λZ). In (4) we made the small angle approximation, neglecting corrections of order (x/Z) 2 and (y/Z) 2 .
For both pixels in our correlation measurement we define coordinates r, φ and z as follows: and a similar set of definitions for pixel 2. Let I(x, y; θ) be the (time integrated) intensity recorded in one measurement at detector pixel coordinates x and y and tilt angle θ. In the experiment we measure the following correlation function: C(x 1 , y 1 ; x 2 , y 2 ; θ) = cov (I(x 1 , y 1 ; θ), I(x 2 , y 2 ; θ)) .
The covariance is estimated in the usual way, from averages of intensities and their products taken from a set of diffraction patterns at fixed tilt. Let I p (x, y, z) be the intensity (squared Fourier magnitude) of the particle density as a function of the scaled wave vector components in the body-fixed coordinate system. Using (4) and the definitions (5), the covariance is expressed in terms of the particle intensity as cov (I(x 1 , y 1 ; θ), I(x 2 , y 2 ; θ)) /N (θ) = (7) where · · · α is an average over the uniformly distributed particle rotation angle α and N (θ) is the number of particles intercepted by the incident beam. In this equation the definition of the particle intensity I p has absorbed several factors, including beam fluence, that we assume are kept constant. In the following we further assume that the particle number density on the substrate is constant so that (again omitting a constant factor) we have N (θ) = 1/ cos θ.
Rotating both the tilted substrate and the detector by π about the beam axis should not change the correlations. This implies the following symmetry property of the correlation function: It is therefore sufficient to only collect data for positive tilt angles.

Data reduction
As a first step in reducing the data we represent the particle density in terms of its angular harmonics with respect to rotations about the body axis parallel to the substrate normal: Using this definition in (7) and the reality property of I p we can evaluate the average over α and obtain We can decouple these equations by averaging both sides of (11) against the factor e −im ′ (φ 1 −φ 2 ) in such a way that φ 1 − φ 2 is uniformly sampled on an interval of 2π radians. Moreover, the decoupling is especially simple if we can do this while keeping the parameters r 1 , z 1 , r 2 and z 2 all constant.
The desired decoupling is achieved by a one-parameter average for which we use either φ 1 or φ 2 , depending on the value of When |A| < 1 we use φ 1 as the parameter; otherwise, by switching the labels 1 and 2 we restore |A| < 1 and the parameter becomes φ 2 . In the following we assume the labels are such that |A| < 1, making φ 1 the averaging parameter. Taking the ratio of equations (5a) and (5c) for both labels we obtain an equation for the second azimuthal angle in terms of the first, where we take the branch that is continuous with φ 2 → 0 in the limit A → 0. In the opposite limit, is monotonic and ranges over exactly 2π radians whenever φ 1 does. Since we perform the average with respect to the parameter φ 1 , we need to include the Jacobian in order that ϕ(φ 1 ) is uniformly sampled.
Equations (5) and the companion equations involving label 2 determine how the detector parameters x 1 , y 1 , x 2 , y 2 and θ need to be varied along with φ 1 and φ 2 (φ 1 ) so that the body-fixed parameters r 1 , z 1 , r 2 and z 2 are kept constant. Comparing (5a) and (5c) (and their companions) we see that the tilt angle must satisfy where compatibility follows from our definition (13) of φ 2 (φ 1 ). We take the branch of θ(φ 1 ) that lies in the range |θ| < π/2, the center of which (θ = 0) corresponds to normal beam incidence on the substrate. Formulas for the other detector parameters are summarized below: Performing the one-parameter average of the correlation function we arrive at the decoupled equations: Before we study equation (18) in greater depth we examine the second branch of the solution to (13) given by Switching to the new parameter we find that the previously defined functions are changed in simple ways. Starting with we find so that (up to an irrelevant 2π shift) the averaging phase for the second branch is the same as for the first branch when expressed in terms of the parameter φ ′ 1 . The same is true of the Jacobian: Finally, using (16) and (17) one easily verifies the tilt and detector parameters for the second branch simply undergo a reversal of sign: The correlation average for the second branch is thus identical to (18a) with all the arguments of the correlation function replaced by their negatives. The second branch therefore provides no additional information once symmetry (9) is taken into account.
In equation (18a) we average the correlations C, between various pairs of pixels on the detector and various tilts, for particular harmonics m and define a correlation C m intrinsic to the particle in that its arguments are spatial frequencies r and z in the body-fixed frame. The two pixels move on circular arcs in the average, since equations (5) and their companions for pixel 2 imply where the right sides are constant in any average. The left sides in these equations represent the magnitudes of the wave vector changes seen in the laboratory while the right sides are the magnitudes of the spatial x y Figure 2: Pairs of detector pixels (labeled 0-12) in one averaging orbit (18) for the case r 1 = 6, z 1 = 5, r 2 = 9 and z 2 = −6. Lines through the circles show the tilt angle θ, where we have used symmetry (9) to keep the range of tilts between 0 and π/2. The tilt is smallest at pairs 3 and 9 and reaches π/2 (grazing incidence) at pairs 0, 6, 12.
frequencies in the particle responsible for those changes. Figure 2 shows the orbits of two pixels and the corresponding tilt in an average for one particular choice of r 1 , z 1 , r 2 and z 2 . Since all orbits visit tilt angle θ = π/2, where the beam is parallel to the substrate, the correlation function needs to be interpolated over the inaccessible range before the average (18) is evaluated. By using the symmetry (9) it is only necessary to acquire data for θ between 0 and π/2.
In the final step of data reduction the angular harmonic amplitudes I m (r, z) are extracted from the averaged correlation function C m . For fixed m, the function C m (r 1 , z 1 ; r 2 , z 2 ) is a Hermitian matrix with row indices r 1 , z 1 and column indices r 2 , z 2 . The content of (18b) is the fact that this matrix has rank one and its unique eigenvector with non-zero eigenvalue, up to an overall phase, is I m (r, z).
The procedure for extracting I m (r, z) is similar to principal component analysis and tests the integrity of the data prior to the phase reconstructions that follow. First a set of radial (r) and longitudinal (z) spatial frequency samples are selected that, by (26), lie within the beamstop and edge of the detector. The density of the samples is set by the speckle scale of the single particle diffraction pattern. Next, the matrices C m are calculated by averaging correlations for the chosen r, z samples according to (18a) and for |m| ≤ M/2, where the maximum harmonic M is set by the angular speckle scale at the edge of the detector. For each matrix C m one then obtains the dominant normalized eigenvector V m (r, z) and corresponding eigenvalue λ m and forms the estimates From the relative magnitudes of the subdominant eigenvalues of C m one can assess the quality of the data at the chosen resolution. If these are significant, then it may be necessary to work with lower resolution blocks of C m , reduce M , or fix the interpolation of the correlation function at θ = π/2 until single-eigenvector dominance is realized for the C m . On the other hand, if no amount of data truncation satisfies this condition then the data are too flawed, e.g. by background, to proceed any further.
We use the ratio of the largest eigenvalue magnitude to the sum of the magnitudes of all the eigenvalues of C m as an internal measure of the consistency of the reduced data, which, as expressed by the notation, is the ratio of the spectral and trace norms of the matrix C m . In the absence of noise and other complications σ m = 1.
Since the phases of the eigenvectors V m in (27) are undetermined, our information about the intensity is incomplete at the level of M/2 unknown phase angles. These angles have to be reconstructed on the basis of additional constraints before the 3D intensity of the particle can be synthesized from its angular harmonics I m (r, z). Since M scales only with the diameter of the particle, and not its volume, this missing information is relatively minor in comparison with the many phases that must be reconstructed to obtain the particle density from the intensity. An algorithm for reconstructing both types of phases is described in section 4.
In contrast to the m = 0 harmonics I m (r, z) that require correlation analysis, the m = 0 harmonic is obtained by a straightforward average of the intensity and is therefore much easier to acquire. There is also no phase ambiguity in I 0 (r, z). Instead of (6) we define A(x, y; θ) = ave (I(x, y; θ)) and (7) is replaced by We can aggregate the three parameter data on the left in a way that tests the two parameter model on the right. One method is to use the earlier parameterization (16) and (17) to define harmonics such that A 0 (r, z) becomes the estimate of I 0 (r, z) and power in the nonzero harmonics indicates data inconsistency.

Data reduction for axial projections
The problem of reconstructing an axial projection from data at tilt angle θ = 0 has been studied extensively by Dilano and co-workers [3,4]. In this section we consider this special case and recover the main result of the previous studies. The final step of our data reduction, however, is different and is the starting point for the improved phasing method described in section 4.1.
By equations (5) the body-frame wave vector z vanishes at tilt angle θ = 0. The corresponding limit of the particle intensity harmonic function I m (r, 0) provides information about the projected particle. Equation (11) reduces to C(x 1 , y 1 ; where the pixel coordinates are now simply and similarly for pixel 2. Specifying the angle of pixel 2 as φ 2 = φ 1 −ϕ, the right side of (32) is independent of φ 1 and we may average the data on the left side with respect to this angle. Decoupling with respect to the harmonic m is now accomplished by a simple integral with respect to ϕ: As in the general case, (34b) implies that the intensity harmonics are given by where V m (r) is, up to a phase, the unique eigenvector of the Hermitian matrix C m (r 1 , 0 ; r 2 , 0) with nonzero eigenvalue λ m . The integrity of the data is assessed using the ratio σ m (28) as in the general case. The harmonic I 0 (r, 0) is just the z = θ = 0 limit of (31) and corresponds, for m = 0, to a simple angular (powder) average of the intensity.
We illustrate data reduction for the axial projection of a simple test particle whose intensity is shown in Figure 3. Uncorrelated, uniform amplitude and normally distributed noise η was added to the true intensity correlations C to simulate the effects of background. Our signal-to-noise ratio in this model is defined by where the root-mean-square amplitude of C is evaluated as a uniform average over all the pairs of measured pixels.
After performing the angular averages (34a) we extract the angular harmonics I m (r, 0) as the dominant eigenvectors of the matrices C m . The result of this for the test particle data with zero noise is shown in Figure 4. That the harmonics I m (r, 0) have a lower triangular support is a direct consequence of the uniform size of speckles in the intensity I p . Because our algorithm for extracting normalized eigenvectors used an unspecified convention for the phase, there is no significance to the global phase of any of the rows of the harmonics plotted in Figure 4 (an arbitrary color-wheel rotation may be applied to each one). The true phases α m have to be determined in the final particle reconstruction stage.
When noise is added to the intensity correlations the harmonics I m (r, 0) shown in Figure 4 become noisy as well. This is well diagnosed quantitatively through the principal eigenvector dominance measure σ m and is shown plotted for our test particle in Figure 5 for three noise amplitudes. We will see that for our particular test particle as few as 10 reliably extracted harmonics (SN = 200) is nearly sufficient for density reconstruction while less than 5 harmonics (SN = 5) is insufficient.

Particle reconstruction
As remarked in connection with equation (27), additional constraints are required to determine the M/2 phases of the intensity angular harmonics I m (r, z). Non-negativity of the intensity is clearly one constraint that can be used. This constraint, however, is weak in comparison with the much stronger constraint that the intensity is the squared Fourier magnitude of a particle density having compact support. And since the particle density is our actual goal, we should approach the reconstruction of the M/2 intensity phases as part of a single process applied to intensity and density jointly. As we will show, the joint reconstruction is a relatively straightforward application of a general procedure once the constraints and the projections to them are set down [7].  Figure 4: Result of the data reduction stage in the reconstruction of the axial projection of a test particle density. The plot shows the complex angular harmonics I m (r, 0) derived from the principal eigenvectors of the angularly averaged correlation matrix (34a), for each m, of the intensity in Figure 3. Magnitude and phase are rendered as brightness and hue, respectively. There is a global phase ambiguity for each m (row), the angles α m , that the data reduction does not determine.
There are three constraints that apply to the pair {ρ, I}, whereρ is the Fourier transform of the particle density ρ and we have dropped the subscript p on the particle intensity I: data : 2π 0 dϕ 2π e −imϕ I(r sin ϕ, r cos ϕ, z) = e iαm I m (r, z), compatibility : The first two apply toρ and I individually and can therefore be combined into a single constraint projection. Projectingρ to the support constraint is no different from how this is done in the standard phase reconstruction problem. Firstρ is Fourier transformed, the resulting density ρ is set to zero outside the support region S and also in its interior when the density is negative, and the result of this is inverse Fourier transformed to give the projectionρ ′ . The data constraint involves the harmonics I m (r, z) given by the data reduction and the unknown phases α m . Details of the corresponding constraint projection, which in effect determines the phases α m , are given below.
With the support and data constraints combined, the compatibility constraint becomes the second constraint in the general scheme where a problem is divided into just two easy constraints. This division of the reconstruction problem is similar to the strategy known as divide and concur [8]. The variablesρ and I are to a large extent redundant representations of the same information and by making them independent we facilitate the projection to the "divided" constraints (support and data). In our case "concurrence" of the redundant variables is imposed by the compatibility constraint, which takes a nonlinear form.
With any projection to a particular constraint set comes the understanding of the distance in the space of variables that the projection minimizes. In our case the variables are of two types: complex Fourier amplitudesρ(x, y, z) and non-negative intensities I(x, y, z), both sampled on the same Cartesian grid. The Euclidean distance in this space of variables, between reconstructions {ρ, I} and {ρ ′ , I ′ }, should be written where the parameter w addresses the fact thatρ and I first of all have different units. Since the compatibility constraint satisfied by the projected variables, at each grid point (x, y, z), is the equation we see that w has the same units as I. Although the precise value of w has a direct effect on the projection {ρ, I} → {ρ ′ , I ′ }, the projected values will always satisfy (39), which is independent of w. The performance of the reconstruction algorithm, on the other hand, can be adversely affected when w is either very small or very large. In those limits one of the variable types,ρ or I, will be much more compliant than the other and as a result the constraint associated with it (support or data) will be explored more rapidly than the constraint imposed on its stiffer counterpart. We obtained good results in the numerical experiments described in section 4.1 when w was set to the maximum measured value of I 0 (r, z), near the minimum of q 2 = r 2 + z 2 , and making w larger or smaller than this by a factor of 2 had little effect on convergence.
If necessary w can be allowed to vary with q while only slightly complicating the projection to the data constraint, the only other projection dependent on w.
Projecting to the compatibility constraint (39) reduces to a one-parameter numerical optimization at each grid point (x, y, z). Given input variablesρ(x, y, z) = ve it and I(x, y, z) = I, the output of the projection is given byρ ′ (x, y, z) = v ′ e it and I ′ (x, y, z) = I ′ where v ′ and I ′ satisfy and minimize The resulting cubic equation for v ′ always has a unique, positive solution which can be tabulated in advance for efficiency.
The nontrivial part of projecting to the data constraint involves only those spatial frequency combinations (r, z) whose magnitude q lies between q min at the beamstop and q max at the edge of the detector. For q < q min , the projection merely copies the input intensities into the output intensities, while for q > q max (in the corners of the intensity grid) the intensities are set to zero.
We uniformly sample the half-annulus defined by q min < q < q max and r > 0 on a Cartesian grid with equal spacings for r and z, and finer than the speckle scale. Every sample (r i , z i ) defines a circle of radius r i in a plane at level z i of the intensity. The first step of the (non-trivial) data projection is extracting the angular harmonics of the input intensity by Fourier transforming the intensity on these circles: where the angular samples ϕ j = 2πj/M are equally spaced and M is chosen to match the angular speckle scale at r = q max . From the data reduction stage we know these harmonics up to an overall phase: All that remains is to determine the angles α m by minimizing The metric of (44) is equivalent to the metric (38) by Parseval's theorem. Since α 0 = 0, the projection simply replaces I 0i by the powder averageĨ 0i . For m = 0 the minimization of (44) yields So far we have defined two projections where P D , the "divided constraint projection", imposes the support constraint (37a) onρ and the data constraint (37b) on I, while the "concurrence projection" P C imposes the compatibility constraint (37c) between ρ and I (while ignoring the support and data constraints). As distance minimizing operations, these are in a sense "greedy" moves to their corresponding constraint sets. Solutions {ρ sol , I sol } to the reconstruction problem have the property that they are fixed by both projections, and in particular,ρ sol is simultaneously consistent with the intensity correlation data and a compact support. There are provably convergent algorithms for finding fixed points in a general projection setting when the constraint sets are convex. Because in our case constraints (37b) and (37c) are non-convex, convergence cannot be proven and we are forced to use algorithms that have a demonstrated record of discovering solutions even for difficult non-convex problems. We use the difference map iteration [7] {ρ, which is equivalent to Fienup's hybrid input-output iteration [9] in the standard reconstruction problem with support and Fourier magnitude as the two constraints (P C → P S , P D → P F ). For our initial iterates we generate a random positive density ρ (uniformly distributed contrast) on the support and from this obtainρ and I = |ρ| 2 , a configuration that only satisfies the compatibility constraint. Iterations are performed until the metric ∆ of the updates (38) fluctuates about a steady state and shows no sign of further decrease.

Reconstruction of axial projections
The specialization of the reconstruction algorithm to axial projections amounts to simply setting z = 0 in the equations of the previous section. In the flat Ewald sphere limit we can also take advantage of Friedel symmetry, thereby reducing by a factor of 2 the work in projecting to the data constraint (the angular Fourier transforms act on a periodic function on the domain 0 ≤ ϕ < π).
The advantage of simultaneously reconstructing intensity and particle density is evident from the frame by frame development of these in the early iterations of the difference map. There is a "locking" of the : Time series of the difference map error metric ∆ in test particle axial projection reconstructions for three values of noise. The corresponding particle densities are shown in Figure 8.
angles α m , that determine the intensity reconstruction, at about the same point in time that the density is well localized in the support, a necessary condition to have speckles in the intensity. In the subsequent refinement iterations the particle density and intensity often executed significant (synchronized) rotation. We deliberately chose a rather loose support to enable this freedom in the reconstruction. Figure 6 shows three frames in a reconstruction with low noise (SN = 10 4 ). The difference map error metric is plotted in Figure 7 for three levels of noise and not surprisingly shows that the fixed-point property (∆ → 0) is strongly compromised when the noise is high. The corresponding reconstructions are shown in Figure 8.
Given the relatively small number of phase angles required for intensity reconstruction, the weaker constraint arising from the support of the intensity Fourier transform -the density autocorrelation -might avoid the complications in the simultaneous intensity/density reconstruction described above. The latter constraint is a weaker constraint on the intensity in that the Fourier transform is only required to have a particular support, and not necessarily the structure of an actual autocorrelation. We tested this simplification  Figure 9: Difference map error metric when reconstructing just the test particle intensity from lownoise data and the autocorrelation support constraint. Convergence is much slower with this weaker prior constraint than for the joint intensity/density reconstructions (Fig. 7). A well reconstructed intensity (not shown) first appears only after 500 iterations. for our test particle by pairing the projection to the data constraint (37b) with a simple autocorrelation support projection in the difference map scheme. This completely avoids manipulating the particle density and the non-linear compatibility constraint (37c). Whereas this simpler approach did succeed in reconstructing the intensity, the progress of the algorithm, shown in Figure 9, was much slower and more sensitive to noise. And since the density is always the final goal of the reconstruction anyway, this two-stage approach does not appear to offer any advantages over the simultaneous intensity/density reconstruction.