1. Introduction

In imaging applications ranging from astronomy to molecular biology, diffraction-based imaging techniques have been demonstrated to be useful. In general, the diffraction from an optically thin object in the far field at low angles can be approximated as the Fourier transform of the scattering power of the original object. Thus, to reconstruct the object, the inverse Fourier transform should be applied to the collected pattern. However, imaging detectors tend to only record light intensity (amplitude squared), losing the complex phase of the underlying wave.

In crystallographic applications, common in imaging of biological particles like proteins, this phase retrieval problem becomes under-constrained, based on the physical data alone. The so-called Bragg peaks representing constructive interference between all unit cell repeats will not, in general, constrain the content of the unit cell to a single representation, i.e. a single phase assignment. By introducing knowledge of possible biological structures in different ways, the problem still becomes tractable [1].

For isolated particles, phase retrieval based on oversampling relative to a support constraint has become the method of choice, i.e. arriving at a set of phases that is consistent with the region surrounding the particle being exactly 0 [2]. To be specific, the critical point is when the number of magnitude observations are twice the number of non-zero pixels allowed by the support [3]. The resulting equation system is non-linear, and solving it in a minimum residual sense results in a non-convex optimization problem. Since it is in general hard to distinguish between local optima and the true global optimum in a non-convex setting, most solution methods do not make claim on providing global optima, but try to balance the quality of the optimum and the exploration of a reasonable portion of the domain. Popular solution methods use variations of alternating projection techniques, such as Hybrid-Input Output (HIO, [4]) and Relaxed Averaged Alternating Reflections (RAAR, [5]). A concise comparison of various iterative schemes and their similarities can be found in [6].

Several of these iterative methods were developed for conditions of relatively high signal, with no or only moderate noise within the recorded field of view. There is an underlying assumption that there exists a set of phases assigned to the exact observed intensities matching the support constraint. Empirically, the methods have been shown to work with moderate levels of noise as well. Still, this is very different from a sparse regime with pixel detectors able to record individual photons, which is now the reality in e.g. FXI (Flash X-ray Imaging). In FXI, individual biological particles are imaged using femtosecond-duration X-ray pulses delivered using X-ray free electron lasers, such as the Linac Coherent Light Source (LCLS [7]). In addition to the sparsity of the signal, some parts of the pattern might be missing due to detector geometry, detector imperfections, or the need to shield the detector from the central undiffracted beam. Despite these limitations, successful reconstructions have been reported for organic and inorganic samples alike [8–11].

In this communication, we investigate a weaker version of the support constraint, based on the known relationship between the Fourier transform of the original wave, and the Fourier transform of the observed quantity, i.e. the intensities. We show the convex character of the optimization problem of estimating maximum-likelihood intensities given a recorded sparse sampling of the diffraction pattern, a support, and a noise model. A wealth of existing methodology for convex optimization can then be used to identify points arbitrarily close to the true global optimum. We name our approach Convex Optimization of Autocorrelation with Constrained Support (COACS). This is implemented as a separate pre-processing step computing the maximum likelihood estimate of the intensities, given the observed data. These new, “healed” non-sparse intensities are then treated as exact in the ensuing phase retrieval process.

COACS shows similarities to methods such as PhaseLift [12,13] and PhaseCut [14], in that the problem is formulated as a relaxed problem compared to the true phase retrieval problem, in terms of a convex or semi-definite optimization problem. This means that a number of existing, efficient, optimization algorithms and software packages can be used. However, we deliberately disentangle the identification of optimum intensities from finding the phases. This, together with great care in the numerical implementation, seems to partially overcome the numerical stability concerns that led the authors in [13] to conclude that not only would PhaseLift be unsuitable for phase retrieval of oversampled particles with known support, but that the problem itself is intractable to numerical solution without adding further constraints.

The global optimum to the full phase-retrieval problem will not be identified with our approach. By pre-processing the intensities, we create a situation where there should exist a global solution to the original phase retrieval problem with a negligible gap between the two constraints, in contrast to what is generally the case for noisy data. This makes it easier to interpret and reduce the errors in phase retrieval. It should also be noted that, with fully exact intensities, 2D phase retrieval in the oversampled case can be reduced to solving polynomials for the locations of exact zeros in the respective 1D transforms [15]. However, this requirement of exact intensities is so strict that the method fails even with quad precision arithmetic on even moderately sized problems. One step in stabilizing the method would require enforcing the solutions to a linear system to be integral −1 or 1, which is not necessarily cheap. For the 1D case, PhaseLift without trace regularization has been noted to be largely equivalent [16]. Still, it is thus conceivable that with near-perfect intensities, a quick and efficient polynomial-based approach might be possible by solving the individual 1D problems. The relevance of our work is then to transition from the full 2D dataset into the near-perfect intensities needed for the decomposition into separate 1D problems.

We demonstrate a proof of concept solution method using a pre-existing software package for compressed sensing, TFOCS (Templates for First-Order Conic Solvers) [17], which was also used in the original presentation of PhaseLift. We evaluate traditional phase retrieval with and without COACS pre-processing based on TFOCS, realizing marked improvements. We also compare our approach to two existing approaches for modified phase retrieval, where an integrated healing step was introduced [18] and with a tuned Gaussian filter on the area outside the support [19].

In the remaining sections, we first describe our optimization problem in greater detail and in relation to the state of the art, followed by a detailed account of implementation concerns, experiments, results, and their evaluation.

2. Phase retrieval with a compact support

We will focus on the case of 2D imaging of 3D objects, since it is straightforward, yet rich enough to illustrate most of our points. The overall methodology generalizes to 3D, like the iterative phase retrieval used as the final step in [20].

Assume that the complex scattering power of the object is P. Assuming ideal conditions (a plane incident wave, a thin object, a far-field regime, and only considering low scattering angles), the scattered wave X meeting a square detector satisfies the following relation:

Since actual photon detectors use discrete pixels covering a limited region in space, we will treat a discretized, finite case, where P, X (matrices in bold) are 2D discretizations of the central part of their continuous, infinite, counterparts. For practical purposes, we can then consider the following version of eqn. 1:

However, these two formulations are not identical. The discrete Fourier transform assumes that the underlying object is periodic. The angles to which the Fourier approximation of diffraction is valid are finite, but the mathematical model of the approximation assumes a process that is infinite in space. Figure 1 illustrates the effects of discretizing this idealized diffraction pattern of our simulated test particle, with small but non-negligible artifacts reaching the edge of the image when transformed back into real space using the discrete Fourier transform. In order to eliminate these artifacts, and their resulting violation of the compact support assumption in phase retrieval methods, we extend (2) with the Hann window W for our square diffraction pattern with side N.

 

Fig. 1 Simulated particle with and without Hann window apodization. a) Original simulated diffraction pattern with high-density sphere. b) Central region of discrete Fourier transform of a), showing the particle. c) The image in b) with limited range to showcase artifacts outside of the object outline. d) The simulated pattern with Hann window applied, resulting in lower intensity away from the center. e) Discrete Fourier transform of pattern after apodization. Slight blur visible. f) The image in e) with limited range. Artifacts found in c) mostly absent.

Download Full Size | PDF

While the need to apodize the signal at the edge might seem trivial, several high-profile reconstructions have in fact been published with appreciable signal all the way to the edge of the detector, offering scant details on how these artifacts are handled. Among these reconstructions are two separate attempts to identify the 3D structure of the bacteriophage PR772 [21,22], where both attempts result in features that could be the result of aliasing, including a conspicuous lack of intensity in the very center of the particle. The problem is aggravated when large regions of missing data exist in the center of the pattern. While the actual artifacts outside of the compact support are of very high frequency, their overall intensity can be reduced by modifying the low frequency components. This will be at the cost of causing ringing within the support.

A related issue is the point where the Fourier pattern reaching a flat detector is no longer a valid approximation of the true diffraction surface. This is only valid at small angles, and at higher angles one needs to adjust for the curvature. It has even been suggested to actively exploit the curvature to reach a 3D reconstruction from a single 2D shot [23]. The critical aspect to note here is that the curvature correction is dependent on the actual geometry of the experiment. For higher angles, the curvature can become relevant. The need for apodization is true at any angle, though. It can only get practically irrelevant if the signal strength is so weak that the overall noise at the edge of the detector dominates over any putative apodization.

If we can assume a photon-counting detector with uniform quantum efficiency r, our observable quantity will be a sampled diffraction pattern in terms of an integer matrix B, where:

Even a strong diffraction pattern will tend to have a characteristic structure of minima, similar to those from the most trivial slit experiment. Thus, even though the pattern as a whole is strong, the Gaussian or identity approximation might not hold for all pixels within the pattern. These regions will give a correspondingly small contribution to the overall signal of the object, but it is still necessary to resolve them to extract any high-resolution features. Most of the contrast in the image of a typical biological particle at non-atomic resolution will be at the border to vacuum, but most aspects of interest to imaging are not in the low-frequency contour. Rather, the interest lies in the finer details of interior and exterior structure. For example, a roadmap regarding LCLS single particle imaging explicitly states a goal of near-atomic resolution for biological particle shots, defined as 3 Å [24].

Reducing even small sources of noise is also relevant when the missing data region is large. A small change in a low-intensity region, when combined with a loose support constraint, can result in large changes within the missing data region.

2.1. Support constraints

Since a single particle is imaged in isolation, we know that the particle has a compact support. We can formalize this by defining an indicator matrix S, requiring that all elements of P outside the support S are zero (P_S^∁ = 0, where the indicator matrix S^∁ is the complement to S). If we simplify to the case where our observations are exact, the following equation system in terms of P (and its Fourier transform X) is obtained:

Due to the squaring of the elements of X, this system is non-linear. Alternating-projection methods are very commonly used for solving such systems. They were pioneered by Fienup [4], and are still popular in practical applications use, as e.g. implemented in the hawk software program and its associated library libspimage [25]. The name comes from the fact that they alternatingly implement some variation on the support constraint, and the intensity (“Fourier”) constraint. The standard implementation is then to assume that the intensities are exact, like in the case above. Both of these constraints can be varied, to better explore the full solution space and accelerate convergence [4,5], or in order to account for measurement errors in the intensities [26]. More recent attempts include using the Alternating Direction Method of Multipliers for easily enforcing more general constraints in both spaces [27]. The exact probability structure of the intensities is rarely used.

The quantity that is sampled to form the directly observable B is Y = X ⊙ X̄. In crystallographic applications, the Fourier transform of Y is referred to as the Patterson function. Using the convolution theorem, we can characterize this as:

2.2. Semi-definite and convex relaxations of phase retrieval

The attempts presented above share the property that they do not claim global convergence. Recent advances apply semi-definite and convex programming approaches to the full phasing problem, finding low-rank factorizations for matrices in order to produce the phases. Such methods include PhaseLift [12] and PhaseCut [14]. These approaches claim to produce unique and globally optimal solutions under favorable conditions with a vanishing probability of not doing so. Example of measurement regimes discussed include observing an object using alternating binary masks, or from unique angles randomly distributed over a sphere [12]. While such conditions are possible to achieve in some situations, they would be challenging to transfer to the single particle X-ray diffraction experiments.

PhaseLift considers the rank-one matrix Z = xx^T, where x is a vector representation of those elements in our diffraction pattern X that are allowed to be non-zero (e.g. a rectangular pixel matrix flattened to a vector), and T indicates the transpose operator. For large supports, in pixels, this will result in very large matrices. Even with an efficient implementation, one expects to do a full 2D FFT per element in x. It is possible to formulate an alternate operator based on the Fourier transform that goes from such representations to the diffraction pattern. The problem can then be solved in terms of satisfying the observed data with a suitable matrix in this space, while also minimizing the rank, i.e. the number of distinct vectors x_i that are needed to represent $\mathbf{X}=\sum x_i{x^{\prime }}_i$. A common way to express the desire of a low-rank factorization of a matrix in convex optimization is to include the matrix trace (the diagonal sum) as a term in the objective function. It is worth noting that with no other terms in the objective function, scaling a positive semi-definite matrix by an arbitrary factor in the range (0, 1) will result in another positive semi-definite matrix, with a lower trace. Thus, the trace term promotes low-rank factorizations as well as a lower overall matrix norm. In [13], this is expressed in their eqn. (2.7) as:

This argument is supported by additional claims about numerical stability. However, in the case reported, noiseless data was supposedly used. Still, the error in the Frobenius norm for the reconstruction, ‖A(X) − y‖₂/‖y‖ was consistently above 0.005. In noiseless conditions, we argue that an error of 0.5% is huge and indicates failed convergence.

Furthermore, it is not clearly specified whether the converged result was a rank 1 factorization or not for these specific experiments. The error metrics presented do not use the converged value of X, but rather a modified matrix where the dominating eigenvector was extracted, renormalized to match the L₂ norm of y, and then recomputed to form the full matrix X (compare section 4.2 in [13]). This normalization hides leakage into secondary eigenvectors, as well as the effect of the trace constraint suppressing overall intensity.

Furthermore, the authors only compared their results to the most trivial traditional alternating projections method, error reduction. In the preparation of our simulations, we did attempt using PhaseLift in TFOCS in a way very similar to what is outlined in [13]. Without changing the constraints, we were unable to find a value of the trace regularization parameter λ which promoted a low-rank factorization, without adversely shifting the overall intensity of the reconstruction. When renormalizing the reconstructions, they were still inferior to reconstructions obtained by HIO, even on noiseless toy size examples. Recently, other authors have also shown that PhaseLift indeed fails to converge in a large number of realistic phase retrieval imaging cases, where alternating-projections methods do produce reasonable results [28]. There are also other semi-definite and convex relaxation methods for the generalized phase retrieval problem. A review of some of these can be found in [29], some of which claim better numerical stability.

Our analysis of the support constraint specifically, will allow an even more favorable numerical treatment as well as a reduced computational load, while still employing a convex framework with its supposed guarantees of a global optimum to the relaxed problem, if numerical problems do not arise.

2.3. Comparison of TFOCS to ER and related methods

For a long time, so-called interior point methods have been the prevailing choice for many families of convex optimization problems [30], with free [31] as well as commercial [32] solvers, and more high-level modeling environments [33]. These methods explicitly form the Hessian of a problem being solved, with a number of elements equivalent to the square of the number of variables. They tend to rely on sparseness in the dependencies between variables. Unfortunately, whereas the support constraint is one type of sparsity, the Fourier operator itself is dense. In some solver implementations, the operator also has to be formed explicitly as a matrix, rather than the far more efficient Fast Fourier Transforms. The benefit of interior point methods is that they tend to converge very quickly, in terms of the number of iterations.

We have chosen to rely on first-order methods, where only the gradient is computed. The TFOCS [17] package is a Matlab toolbox, and is presented as a set of templates from which one can build tailored solvers. There are two main forms for formulating problems in TFOCS, and we have used the one that corresponds to the eponymous function tfocs:

The linear operator 𝒜 does not need to be implemented as a matrix, but can be expressed in code, assuming it satisfies certain specific criteria. We can thus rely on an actual efficient Fourier transform implementation. The function f needs to be smooth and differentiable, with an explicit implementation of the gradient. For h on the other hand, it is enough that it is possible to evaluate the proximity operator, i.e. that it is possible to (quickly) solve:

This is sometimes referred to as the implementation of the function being prox-capable. For a projection, the proximity operator and the projection are identical. In this context, the projection can be seen as an indicator function that takes the value of 0 for admissible points and inf everywhere else.

It can be helpful to note that the traditional error reduction problem can also be expressed in the form of (12). In this case, though, both functions are projections. Furthermore, the practical implementation consists of alternatingly fully implementing one constraint or the other. TFOCS, on the other hand, will make iterative updates based on the gradient of f and the proximity operator of h, with various specific algorithms. Beyond the most trivial gradient descent algorithm, all optimization algorithms implemented in TFOCS are capable of converging faster than a trivial alternating projections approach, since the search direction applied is modified to take into account the history of previous iterations, whereas the two projections in ER or HIO can almost cancel each other out, with only a small net step per iteration. The demo_alternatingProjections.m script in the TFOCS source distribution explicitly illustrates this point. Similar accelerated approaches are used in e.g. stochastic gradient descent training of deep learning networks, although their behavior is not proven outside of a convex context. The β coefficient in HIO provides a mean of acceleration, but it is far more static and crude.

3. Existing methods for intensity healing

The concept of intensity healing is not new in itself. Previous approaches have tended to focus on either the missing data problem, in terms of e.g. a beamstop, or the problem that the recorded intensities are noisy.

For the first problem, several methods focus on the presence of so-called missing modes, that is, linear combinations of pixel values that are only to a negligible extent affected by the Fourier constraint as well as the support constraint. The most obvious such mode is a Gaussian function, since the Fourier transform of a centered Gaussian is another centered Gaussian, with inverted standard deviation. If the size of the support and the size of the missing data region allow a Gaussian with negligible tails to be fitted in both spaces, this mode can increase to very high magnitudes. More sophisticated modes can also be computed. Foundational theory for determining missing modes in the discrete case can be found in [34]. More recent practical applications in the context of phase retrieval include [35] and [10].

We do not want to downplay the overall importance of missing modes, but we note that they are less relevant if the size of the central speckle of the transform of the support is larger than the size of the missing data region. In this case, the inverse relationship between the size of features in the two spaces mean that no self-contained mode can fit unaffected in both spaces. Thus, if the central speckle of the object is visible, the missing mode problem is reduced assuming an appropriately tight support is defined. This is demonstrated in practice in our simulation results.

Inferring the correct intensities in a missing data region that is wide, but not covering all of the central speckle, can still pose a practical challenge. High levels of noise in the Fourier signal can also cause unwarranted artifacts in the phase retrieval process. One approach that tries to handle these issues is the iterative algorithm Oversampling Smoothness (OSS) [19]. This is presented as a modification to HIO, where a Gaussian apodization is applied in Fourier space, but the Fourier transform of this apodized version of the intensities is only used for those pixels that are outside the support in real space (eqs 2–4 in [19]). Furthermore, the bandwidth of the filter is adjusted gradually to penalize high-frequency oscillations outside the support more strongly in later iterations. According to the authors, this filter reduces phasing process oscillations caused by noise at higher frequencies.

There have been previous attempts to use autocorrelation support constraints, similar to what we formulated in (10). One of the most successful examples of this is found in [18], where the approach is termed “speckle healing”. Speckle healing is performed as an outer iteration process in phase retrieval, performed every 1,000 iterations. At this point, the current intensity estimates for the missing data region are extracted and combined with the observed intensities. The projection defined by (10) is applied once, resulting in new intensity values everywhere. For the non-missing region, these new intensity values are then treated as the exact truth during the ordinary phase retrieval for another 1,000 iterations.

The speckle healing approach is attractive, in that the speckle healing and phase retrieval processes are connected. The approach ensures that the phase retrieval algorithm does not encounter a sparse or noisy pattern at high frequencies. Rather than partially ignoring the data at high frequencies, like in OSS, the data is transformed into a healed form which should be closer to satisfying the constraints of the phase retrieval process.

However, a single projection application of (10) fails to account for the fact that all intensities should be non-negative. The auto-correlation constraint is also not directly used to infer the intensities in the missing data region. The exact shape of the speckles in the observed region will be constrained by the current intensity estimate from the phase retrieval process in the missing data region. Convergence in the missing data region can be very slow if several such outer iterations are needed in order to make the missing and observed region match properly.

4. Intensity healing using COACS

Based on the structure of the probability mass function in (5), (10) is a likelihood-optimization problem with a linear transform into independent variables. The previous quadratic term for X has vanished, at the possible expense of the self-convoluted version of S, Ŝ, imposing a less strict constraint on the solution.

If one would assume a Gaussian noise structure, the resulting problem from (10) is an ordinary least squares problem, which is itself a convex problem, since the quadratic residual is convex. A solution to this system can be determined by solving the over-determined homogeneous linear equation system based on the Fourier operator for those pixels that are in the complement of the self-convoluted support. This is what was done for the “speckle healing” in [18], since the direct application of the projections amounts to minimizing the least-squares error in the other space. Similar observations of solving over-determined linear equation systems have also been done by others, e.g. in [36].

For a more accurate solution, one will however also need to add the constraint of all intensities being non-negative, which necessitates specialized convex optimization solvers, rather than the most basic least squares routines. Negative intensities are impossible to convert to the scalar amplitudes needed for the original phasing problem. Truncating negative values to zero, as well as taking the absolute value, will cause artifacts. Therefore, we use the TFOCS model equation (12), with the likelihood model from (5) being expressed as f, the support constraint in (10) as h, and 𝒜 representing the 2D Fourier transform. These are implemented as custom functions, for numerical accuracy reasons outlined in section 4.1.

4.1. Numerical concerns

As noted by the authors in [13], numerical stability and accuracy concerns are of great importance in the phase retrieval application, especially when relying solely on the oversampling of a single pattern. Although the point of doing a separate step of intensity healing is to make the phase-retrieval problem more well-posed, these numerical concerns are of importance for intensity healing as well. The intensity in a diffraction pattern of, say, a uniform sphere, will decay in a manner roughly proportional to q⁴, where q is the scattering angle [37]. The dynamic range needed for a problem in intensity space is higher than in the amplitude space, due to the quadratic relationship between amplitudes and intensities. Since the Fourier transform adds terms for all frequencies for each resulting pixel, and terms are of varying sign, loss of precision due to terms of opposite sign canceling out can be paramount. In any iterative approach, the numerical accuracy needs not only to hold for the final result, but it is also necessary for properly evaluating the change undertaken between iterates.

The steps outlined below try to consistently reduce the chance of:

The third cause interplay with the first two. If the problem is formulated in such a way that values in some domain are large, taking a very small step can cause both truncation errors and loss of precision, to the point where the desired step change between two iterations is rounded to exactly 0. This does not mean that convergence has been reached. With an adaptive step size, later steps, away from the domain boundary, can be large.

The specific places where this has been taken into account can be summarized as:

It should be noted that without these steps, one arrives at a seemingly mathematically equivalent problem formulation that still results in erratic numerical behavior in practice. The specific nature of these steps also relate to design choices in the existing TFOCS solver, including the fact that differences between successive iterates are computed, and that the objective function is handled as a scalar, rather than considering the individual per-pixel contributions. On the other hand, these properties are in no way unique to TFOCS.

4.1.1. Lipschitz backtracking structure

TFOCS can be described as an alternating-projections method tracing modified gradients, giving far accelerated convergence compared to a fixed-step gradient descent. For more details on the overall design of the package, we refer to [17, 38]. One crucial aspect of the algorithm is a backtracking approach to identify the proper step size. The backtracking is common to the various specific solver schemes implemented in TFOCS. The methods rely on the gradient for f being Lipschitz continuous, i.e. that the following relationship is satisfied for any vector d:

The global bound on the Lipschitz constant L can be very restrictive, resulting in small step sizes. The convergence of the methods only relies on the rate of variation in the local vicinity currently visited. Therefore, TFOCS uses an adaptive estimate of a local Lipschitz constant for the gradient of f.

As elaborated upon in [17], there are multiple ways to estimate a bound on the Lipschitz constant based on locally evaluated gradients and function values. Since the problem is inherently related to evaluating the difference in function value between two points, there is a great risk for loss of precision. The original TFOCS approach is therefore to use two bounds, one of which will give less conservative estimates, while also being more sensitive to the details of floating-point math, and another bound that is generally more conservative, but numerically safer.

The TFOCS code automatically detects whether the relative difference in function value between the two points is small, and then switches to the more conservative bound. The conservative method will then be used until the next “restart” of the TFOCS algorithm.

We have made two modifications to this logic. The first modification is to move the evaluation point for choosing which bound formulation to use. The original TFOCS code performs one step and then, retrospectively, checks whether the difference was small. With this design, a single step in the iteration process might be taken using an accuracy-deficient Lipschitz estimate. This turned out to be enough to cause drastic divergence for us in some scenarios. We instead put this test earlier. We also allow the iteration process to leave the conservative regime, without a restart. Furthermore, we consider not only the difference in function values |f(x) − f(y)|, but also the difference |x − y|. If the step itself is small relative to the magnitude of the values, truncation errors as well as loss of precision are possible.

The total change is less than 20 lines of code in a single function in the package. It has been submitted to the authors as a suggested change. This, together with using the “no-regress” restart feature gives adequate optimization performance without divergence or singularity issues.

Despite these changes, we have concluded that a fixed step size, based on the known Lipschitz bound for f, can work well for high accuracy. Since the true solution will in general be close to 0 for many pixels, which is the region with the highest Lipschitz constant, a local bound for L will tend to be close to the global one anyway.

4.2. Expressing the solution variables as a perturbation

The Fourier transform is a linear operator. Thus, the relationship between Y and P̂ can be expressed in terms of a perturbation of the relationship between Y⁰ and P⁰ :

To fully realize the benefits of the perturbation formulation, the corresponding shift needs to be done in the computation of the f term of the objective function as well. First of all, consider the negative Poisson log-likelihood parametrized on the rate parameter, with the observed intensity k_i for pixel i. Since it is a likelihood, the factorial denominator in the Poisson probability mass function can be safely ignored.

This will then easily sum to a large value, resulting in the issues we want to avoid. What we can do to improve this is to use Y⁰ as a reference point in both domains, and only compute the difference.

Note that just computing f(Y⁰ +Y^⋆)− f (Y⁰) does result in the same value from a mathematical standpoint, but its numerical implementation results in a serious loss of precision due to computing the difference between two large sums. Rather than computing the difference of the sums, we compute the sum of the individual differences.

This function is a perturbed version of the smooth_logLPoisson function that is already included in TFOCS. To further accelerate convergence, especially for pixels with zero observed photons (where the log-Poisson optimum lies exactly at 0, the border of the non-negativity constraint), we replace the log-Poisson with a linear extrapolation below a border value l, and add a quadratic function scaled to still keep the minimum at 0 for the zero-photon case. Hence, negative values are allowed, but strongly penalized. This modified version of f can be expressed as:

Naturally, the second case can be simplified, but these three terms form the motivation. We evaluate the original function at the border of our barrier, add a linear extrapolation to keep the gradient smooth, and a quadratic function that is zero at the barrier and scaled in such a way that the quadratic and the linear extrapolation cancel out at Bi.

Finally, in our implementation, we have chosen to implement the apodization factor inside our final version of f, rather than in 𝒜. Since the reason to introduce l is to ensure a smooth behavior amenable to optimization, it makes sense that the effective vector Laplacian (gradient of the gradient) should be of equal magnitude everywhere. The smoothness of the gradient of f is what limits the step size in TFOCS. Therefore, the actual value used for the barrier width is inverted by the window squared, resulting in a more lenient barrier in the high-frequency areas, thus also requiring an overall lower value of l to guarantee high accuracy at all diffraction angles. Likewise, we adjust l_i to account for the vector Laplacian of the Poisson log-likelihood, since that will grow very high close to zero when the number of observed Y_i is non-zero. We thus choose a point l_i for each pixel which bounds the maximum gradient of the gradient of f. Thus, the final form becomes:

W with elements w_i indicates the window used. The window we use is based on (3), but with an extra ε term to avoid singularities in the inversion. Thus, it is not a perfect Hann window, but slightly raised. The proper choice of ε is a tradeoff where additional numerical accuracy is needed to realize the theoretical benefits of a more perfect window. In practice, we have used 0.00025.

4.3. Relaxed non-negativity barrier

The non-negativity barrier width l introduced in the previous section determines the sharpness of the non-negativity constraint. In our implementation, we initiate l to a high value (l = 4). Lower values of l result in a more narrow quadratic penalty. As l → 0, the barrier turns into a strict non-negativity constraint, increasing accuracy, but slowing convergence by reducing the allowable step size when optimization variables are close to 0. In our implementation, the value of l is then decreased by a factor of .5 repeatedly until reaching 2⁻⁶⁶ ≈ 10⁻²⁰. Due to the Hann window scaling of l_i, this means an effective barrier of 10⁻²⁰/ε² ≈ 10⁻¹³ at the edge of the pattern.

This scheme ensures a very quick convergence to the relaxed problem, followed by successive sharpening of the barrier in a continuation scheme, starting off from the previous end result. Let X_i denote the result from iteration i (starting at 0, with l = 2²⁻ⁱ). Since the major trend influencing the results should be the reduced fluctuations around 0, there is reason to believe that a first approximation of the optimum when changing l can be expressed on the form X_i+1 ≈ X_i + d(X_i − X_i−1) when i > 2, i.e. a further movement along the direction of change between the previous two iterations. This acceleration scheme is similar to the ones more thoroughly investigated for the SCD model in [17], but we employ a line search in this direction to find the value of d that results in an optimal decrease in the objective function.

The actual X values used are expressed as perturbations for the accuracy reasons already discussed. Each continuation step with a new l value induces a rapid shift in the optimum, even after the acceleration adaptation. Therefore, a second set of continuations is applied after a fixed number of individual inner iterations of the optimization scheme. The termination of the algorithm and continuation restart allows a new perturbation center to be used, and also the application of another line search acceleration. The inner continuation loop is terminated when the change in the objective function falls below a predefined threshold. The number of iterations is increased for each inner continuation step, since sometimes the TFOCS algorithms require a high number of iterations to find a good search direction, balancing the gradient from f and the proximity operator from h.

4.4. Apodized support term

The support constraint can be simply implemented as a projection to 0 outside of the support. However, this might unnecessarily constrain the model in terms of inducing very small steps. Rather than implementing h as a strict projection, we therefore instead use a proximal operator representing a quadratic penalty on values outside the support with strength $c=\frac{5*{10}^7}{l}$:

This can be extended into perturbation form of the individual terms of the sum, analogous to what was showed for f in previous sections.

However, this form of the support term fails to take into account the fact that we have assumed apodization in our solution. In practical experiments, we could observe significant high-frequency content arising at the border of the support. Eventually, these high-frequency signals vanish, as the model converges, but the process is slow. Since the apodization in Fourier space corresponds to a slight blur in the image space, the same holds for the autocorrelation space. With a blurred autocorrelation, the transition between the zero signal support and the non-zero support should also be gradual.

Therefore, we suggest applying the window W to Ŝ as well.

In practice, similar results were achievable without this addition, but the time usage was much higher. The windowed version will only accelerate convergence by avoiding spectral leakage.

5. Data and experiments

Simulated diffraction data were pretreated with our proof-of-concept implementation of COACS based on TFOCS, called jackdaw (available at https://github.com/cnettel/jackdaw).

Simulations were done using the Condor [39] package. The test particle chosen in the simulations consisted of an icosahedron, with a single very high-density sphere placed on a point off-center on one of the vertices. The shortest diameter between vertices was 20 pixels, with the diffraction pattern recorded on a virtual 256×256 detector. The high-density sphere had a diameter of 4 pixels, with a much higher scattering power per volume than the icosahedron. The presence of a concentrated, markedly non-symmetric feature makes it possible to judge the presence of artifacts in resulting reconstructions, as well as whether the feature is resolved at all. At the same time, a rough approximation of our simulated sample is a sphere, which is appropriate for the types of biological particles one would want to image.

A total of 50 Poisson samplings were created of this particle. A central beamstop was simulated by making a 25×25 square stencil in the center missing. This covered the majority of the central speckle for the icosahedron, resulting in on average 10,000 photons in the non-masked pattern. Intensity-healed patterns were computed using our COACS implementation for this pattern, using a trivial auto-correlation support consisting of a square with a side of 61 pixels. While a Shrinkwrap approach akin to what is common in phase retrieval [40] would be possible, this choice was made to demonstrate the lack of specific assumptions on the structure of the sample, and how it would be possible to process a large number of shots of varying particles with little to no manual tuning.

Phasing was performed using the Python interface to libspimage with the GPU-based HIO phasing method, with the relaxation parameter β = 0.9, identically for COACSed (healed) and un-processed diffraction patterns. The support used was a square with a side of 31 pixels, i.e. slightly larger than the particle itself. 50,000 iterations of HIO were performed, followed by 10,000 iterations of pure ER to find a suitable local minimum. For each pattern, phasing was repeated 100 times, in order to control for the variations due to random initial data. The eventual reconstruction chosen for each pattern was the average of the 10 best out of these phasing results, chosen on the basis of the real space error, i.e. the norm of the remaining signal outside of the support.

A modified version of our phasing approach inspired by the speckle healing (SH) approach was also tested [18]. This was implemented within the libspimage code, but every 1,000 iterations, the observed intensities were replaced by a current guess based on the non-apodized autocorrelation support constraint. Since speckle healing is supposed to implement the same constraint as COACS, this was only attempted on non-preprocessed data.

The additional filtering suggested by the OSS method is supposed to handle the noise due to sparsity at high frequencies [19]. Since OSS is supposed to make a transition between HIO-like and ER-like behavior, we executed it with standard settings for 60,000 iterations with 100 replicates, based on the demonstration Matlab source code published by the authors. We tested this with and without COACS. Finally, we tested the relevance of apodization by running COACS and HIO phasing without the Hann window applied.

6. Results

A qualitative comparison of the results possible using COACS, speckle healing and OSS are given in Fig. 2, for the first 10 of the 50 patterns simulated. For each unique random pattern with approximately 10,000 photons outside the beamstop, HIO with COACS, OSS with COACS, HIO without COACS, OSS without COACS, and Speckle Healing (SH) are shown, as well as two cases without apodization. The high-density spherical feature on the edge is clearly visible in all COACS-healed reconstructions, and the contours of the icosahedron projection are reasonably traced. Even with the additional processing introduced by OSS and SH, the non-COACSed versions are not able to resolve this non-symmetrical feature. Furthermore, the edges of the icosahedron are less clearly defined and the noise outside it stronger. Without apodization, the correct features are less readily identified in the COACS case. While SH is supposed to implement a constraint that is equivalent to the one in COACS, the end result is better than pure HIO, but far from HIO + COACS in terms of overall visual quality.

 

Fig. 2 10 phased reconstructions of sparse patterns based on the simulated particle, each with approximately 10,000 photons, with different phase retrieval schemes compared against the true noiseless reference image. The methods include two cases using our COACS pre-processing, with Hybrid Input-Output and Oversampling Smoothness. The same methods were tried without COACS as well, plus Speckle Healing, since the additional constraint in that method is in theory equivalent to COACS. The final two cases verify the behavior when the apodizing Hann window is not applied. Example of an actual input pattern is shown in Fig. 3. Each picture is an aligned average of the 10 best individual phasings out of 100 replicates for each combination of phasing method and sampled image.

Download Full Size | PDF

The diffraction pattern with various stages of processing is shown for the first random instance in Figure 3. The reconstructed patterns look similar to the original in Fig. 1. The phased COACS patterns and especially the non-phased COACS patterns are clearly superior to the patterns achieved using traditional phasing. There are no perceptible ringing effects from the highly regular rectangular support used in the COACS patterns. Finally, a pattern retrieved using COACS without the apodization is shown. The overall structure is similar, but the fringe locations have shifted, indicating a consistent bias.

 

Fig. 3 Diffraction pattern results for the first random instance (rescaled to undo Hann window apodization). a) Original simulated pattern. b) Poisson sampling with simulated 25x25 beamstop, giving only 10,000 non-masked photons. c) Average diffraction pattern over 10 best phasings. Artifacts at the edges, reproduction of original features dropping off quickly. d) Average diffraction pattern over 10 best phasings based on COACS healing. Features reproduced to higher angles, but still problems at the edge, missing the outer ring. e) Diffraction pattern resulting directly from the COACS healing. Artifacts present in d) are absent. Outer ring corresponding to the one found in original pattern. f) Average based on the separate healing results of all 50 sampled particles. Errors cancel out, resulting in a pattern very similar to a).

Download Full Size | PDF

The crystallographic R factor for each reconstruction was computed relative to the simulated data in the same way as presented in [41]. This is equivalent to the L₁ norm of the difference between the recovered Fourier space amplitudes and the true simulated amplitudes, normalized by the norm of the true amplitudes. R factors together with real-space mean-square errors of aligned reconstructed images, where available, are presented in Table 1. The R factor with no healing is more than 100% higher than what is obtained when phasing with a healing pre-processing step. However, even this result has an error significantly larger than what is obtained from the healing step itself. Thus, the phasing process is still a significant contributor to errors. Since these calculations go all the way to the edge of the detector, where the sampled pattern was exceedingly sparse, the values are relatively high. The MSE results also show marked improvements. R factors and MSE produce results that are consistent with the qualitative visual judgment already presented, with COACS + HIO coming out on top, closely followed by the other COACS variations.

Table 1. Pattern R factors and image mean-square errors.

View Table

For a smaller selection of the methods, we also present the R factor calculated at various radii (in pixels) in Fig. 4. This radial R factor has sometimes been referred to as the R Factor Transfer Function (RFTF) [41]. The non-COACS methods are not able to correctly recover the intensity in the central missing data region (being a square of side 25, this is seen most up to a radius of 12.5). In the overall minimum for the true signal at 90 pixels, due to the shape of the high-density spherical feature in the image, the direct application of COACS gives inferior results due to the more specklish nature of that pattern. Although the relative error is higher at this point, the absolute error is small nonetheless, since the true signal is close to 0.

 

Fig. 4 R factors (normalized relative L1 error) for various radius shells in pixels. Curves are averages over the individually computed results for all 50 simulated particles. Comparison between results based on average phasing of 10 best reconstructions of original pattern, average phasing of 10 best reconstructions of COACS-healed pattern, and using the structural factors from the COACS pattern directly. COACS-healing reduces phasing errors, but the phasing step is still a significant contributor to errors. OSS is competitive with HIO overall, but fails at reproducing correct signal at the lowest frequencies, which also contributes to the higher MSE error for that method. Speckle Healing, which theoretically implements the same autocorrelation constraint as COACS, improves results over pure HIO, but does not come near HIO + COACS. COACS peak at around 90 pixels is due to the R factor being a relative error metric. This is the location of a minimum due to the shape of the small spherical feature. Hence, absolute errors of the same magnitude are amplified.

Download Full Size | PDF

The healing process required approximately 120 minutes per pattern on a 12-core server with Sandy Bridge Intel Xeon cores. In the Discussion section, we outline how this can be improved.

7. Discussion

We have demonstrated superior visual quality, as well as a clear reduction in the R factor and mean-square error when applying COACS as a pre-processing step before phasing. The time usage for the intensity healing is much longer compared to phasing with HIO, 120 minutes vs. less than 10 seconds (for a single reconstruction within the batch of 100). Nonetheless, both perform roughly the same number of iterations (on the order of 10⁵).

It should be noted that the current COACS implementation is a straightforward one in Matlab using the TFOCS toolbox, while the phasing implementation libspimage [25] is efficient GPU-based CUDA code. When using the TFOCS mode for counting the number of operations, it is clear that the number of function and transform evaluations scales linearly to the number of iterations. A preliminary analysis of the computational workload indicates that the Fourier transforms should be the dominating part of the computational load. Thus, an efficient implementation of the TFOCS algorithm should be quite similar to the existing phasing, in terms of work per iteration.

The time difference is due to the difference between a CPU-based implementation in Matlab for intensity healing, and an efficient C++ implementation on GPUs using CUDA. Another avenue would be to use a library such as cvxflow, a high-level tool for expressing convex problems and then using a solver that executes efficiently on the GPU based on the TensorFlow framework. Note that this would not be an actual application of deep learning methodologies to the phase retrieval problem, such as the one in [42], but a matter of employing the effective implementation of computational graphs in TensorFlow for any task heavy in linear algebra, in this case the first-order methods of TFOCS. The OSS implementation we used, based on the sample code in [19], was also a straightforward Matlab implementation. A single phasing using this method consumed 20 minutes. A set of 100 reconstructions therefore uses over 30 hours without additional tuning or parallelization, which makes the 2 hours used by the current COACS implementation acceptable in comparison.

One can also easily improve performance by a factor of four or more with a limited increase in the R factor, by reducing the convergence thresholds within TFOCS. However, this seems to be at the cost of a higher variance in the results. When tested, a few outliers of our 50 random cases would display error metrics distinctly higher, indicating cases which did not fully converge. The phasing process following the pre-processing should supposedly be possible to speed up significantly. In our simulations, 256×256 patterns were used with a 31×31 support. A downsampling to e.g. 64×64 would be straightforward, still provide the sufficient oversampling, and supposedly reduce the work per phasing iteration by more than a factor of 16. We chose not to test this in our experiments to keep the configuration of the phasing algorithm as similar as possible for the cases with and without COACS pre-processing. Without the pre-processing, similar downsampling is not possible due to the missing data region.

We believe that slight additional relaxations like the ones already outlined, and an efficient GPU-based implementation, should make intensity healing achievable at time scales similar to phasing. More importantly, we have demonstrated that a proper numerical treatment allows proper intensity healing and phasing without a carefully chosen tight support or Shrinkwrap procedure. This development holds the promise of allowing routine phasing of recorded experimental diffraction patterns, without tedious manual tuning. With equally tight support constraints, it is also foreseeable that this method will allow better handling of patterns with missing data in some parts of the pattern due to e.g. saturation, the presence of a beamstop to protect the detector from the much stronger non-diffracted beam, or other aspects of detector or experiment geometry. For analyzing real-world data, one will probably want to augment our current Poisson probability distribution with a Gaussian component for values close to 0, in order to better account for electronic noise in imaging detectors.

Another observation based on our evaluation is that the phasing schemes are still the weakest step in the reconstruction procedure. Even the smooth COACS patterns result in the object in the real-space image showing a clear drift during the reconstruction processes, within the support. This is due to a non-ideal phase ramp being induced by the combination of the two constraints, resulting in movement in real space between iterations. The issue of drift has been described multiple times, including in [35], where the authors also several times note the benefit of a tight support to avoid such issues. We did not choose to tighten the support in our experiments, since such tightening almost always includes some level of manual intervention – what is the proper threshold between the outer contour of the object and the surrounding environment? Shrinkwrap approaches tend to work well in practice, but with some level of manual/dataset-specific tuning, especially for data this sparse.

One aspect in how the drift induces artifacts is that the object will not only move within the support, but also reach its border. This results in new artifacts being introduced as the real space constraint implementation tries to remove the signal (especially more advanced implementations than pure error reduction). While the scope of this publication is not to improve phasing per se, new insights can be gained by clearly separating effects due to sparse data, the difference between the continuous and discrete Fourier transforms, and the phasing methodology itself. We also note that a few individual reconstructions look far sharper than the averages shown, although the error metrics do not reliably separate sharp results from less sharp ones. It is our opinion that a proper phasing method should be able to produce R factor values similar to those we report for the COACS method alone in Table 1 and Fig. 4.

Both our preliminary tests and the results of others indicate that PhaseLift can display problems in providing a successful reconstruction in cases where other phasing methods succeed. The resulting matrix can be far from rank 1, which is a requirement for a successful reconstruction. Even when a low-rank matrix is attained, this can be at the cost of a high value of λ, which suppresses the overall intensity of the pattern during the optimization process. Renormalization after completing the optimization process, as was done in [13], reduces the error reported, but it does not cause proper convergence to the correct intensities.

We also note that if PhaseLift is started with a zero matrix, all resulting eigenvectors will be symmetric, when the measurement matrix is the complete Fourier operator. Without an initial asymmetry, the gradient computed from the observations will never result in an asymmetry being introduced. Thus, that method is sensitive to the starting point. It is reasonable that measurement operators that lack the specific symmetry properties of the Fourier operator alone do not trigger this problem in PhaseLift, which might help explain the authors’ far more favorable results in other regimes. For some other problem settings, it seems that setting the regularization parameter λ to zero yields a low-rank factorization, but that is not in general true for oversampled single particles.

Thus, a straightforward application of PhaseLift to an ideal version of the oversampled compact support case fails to converge to a meaningful solution within its own framework, namely one where the matrix X produced by the method is of almost rank 1 while at the same time reproducing the intensities even in noiseless simulated measurements. The conclusions of the authors of [13], that these deficiencies are due to inherent problems in the oversampled phase retrieval problem, are unfounded. There are several other modern schemes for phase retrieval using convex or semi-definite relaxations, some of which claim better numeric performance or tolerance to noise [29]. However, we believe that the specific separation of the intensity problem and the phase problem made by COACS can be a useful addition to the set of available methods, and one that could be readily employed together with existing phase retrieval for actual data from coherent diffractive imaging experiments.

We also want to point out that while apodization proved critical for high-quality results in our method, the question of whether to apply apodization or not is relevant for all phase retrieval for single particles. It is also relevant for periodic structures where the period is not consistent with the field of view of the detector. The problem is only highlighted because of our treatment of all observed intensities as non-exact, with a maximum-likelihood objective, rather than a projection onto the intensities. The OSS approach does apply a degree of apodization, but this is motivated by the need to reduce oscillations in the reconstruction. The varying bandwidth of Gaussians used and the fact that they are only applied to the area outside of the support, do not provide an immediate spatial interpretation of the result equivalent to the 1-pixel blur of a Hann window.

The straightforward structure and generality of the convex optimization formalism in general, as well as the TFOCS library in particular also make it feasible to add additional constraints or regularizations to COACS. Such constraints might include a total variation norm in the autocorrelation or Fourier space to regularize the problem. Especially in exceedingly sparse cases, and with challenging experiment geometries, such additional inspiration from the compressed sensing literature might prove worthwhile. Specific removal or regularization of the missing modes could also easily be added as additional objective function terms. One observation to make about the behavior of TFOCS is also that while the formal requirement on smoothness of f is that the gradient is Lipschitz continuous, the actual backtracking schemes making local estimates of the Lipschitz constant L can display problems if the vector Laplacian is not also smooth. This led to us using a constant step size and a carefully tuned formulation of f even after making fixes to the backtracking logic in TFOCS.

8. Conclusion

We have presented the COACS approach to correct the sampled diffraction pattern based on a support constraint. This approach allows for higher-resolution reconstructions compared to the conventional Hybrid Input Output (HIO) approach, as well as Oversampled Smoothness (OSS) and Speckle Healing (SH). Our preprocessing approach made it possible for us to phase simulated data with a wide support and no specific tuning of the parameters with a high level of detail. We have also discussed why the PhaseLift method fails to converge for the oversampled compact support phase retrieval problem, while other methods succeed. We have also identified several possible future developments, most pressingly to implement a GPU-based version of our approach in order to present more competitive performance in terms of computation time.