Using Convex Optimization of Autocorrelation with Constrained Support and Windowing for Improved Phase Retrieval Accuracy

In imaging modalities recording diffraction data, the original image can be reconstructed assuming known phases. When phases are unknown, oversampling and a constraint on the support region in the original object can be used to solve a non-convex optimization problem. Such schemes are ill-suited to find the optimum solution for sparse data, since the recorded image does not correspond exactly to the original wave function. We construct a convex optimization problem using a relaxed support constraint and a maximum-likelihood treatment of the recorded data as a sample from the underlying wave function. We also stress the need to use relevant windowing techniques to account for the sampled pattern being finite. On simulated data, we demonstrate the benefits of our approach in terms of visual quality and an improvement in the crystallographic R-factor from .4 to .1 for highly noisy data.


Introduction
In imaging applications ranging from astronomy to molecular biology, diffractionbased 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 a projection of the original object.Thus, to reconstruct this projection, the inverse Fourier transform should be applied to the collected image.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.The resulting equation system is non-linear, and solving it in a minimum residual sense results in a nonconvex optimization problem.Thus, most solution methods do not make claim on providing global optima, but try to balance finding optima and exploring a reasonable portion of the domain.Popular solution methods use variations of alternating projection techniques, such as Hybrid-Input Output (HIO, [2]) and Relaxed Averaged Alternating Reflections (RAAR, [3]).A concise comparison of various iterative schemes and their similarities can be found in [4].
Several of these iterative methods were developed for conditions of relatively high signal, with no or only moderate noise.This is very different from a sparse imaging regime with pixel detectors able to record individual photons, which is now the reality in e.g.FXI (Flash X-Ray Imaging), where individual biological particles are imaged in femtosecond-duration X-ray pulses delivered using X-ray free electron lasers, such as the Linac Coherent Light Source (LCLS [5]).In addition, some parts of the image might be missing due to detector geometry, detector malfunction, or the geometry of the experiment itself.Despite these limitations, successful reconstructions have been reported for organic and inorganic samples alike [6,7,8,9].
In this communication, we introduce a relaxed 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, 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.Convex optimization problems lack separate local optima and allow a wealth of existing methodology to be used to identify points arbitrarily close to the true global optimum.We name this method Convex Optimization of Autocorrelation with Constrained Support (COACS).
COACS has similarities to methods such as PhaseLift [10,11] and PhaseCut [12].However, we deliberately disentangle finding 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 [11] 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.
We demonstrate a proof of concept solution method using a pre-existing software package for compressed sensing, TFOCS (Templates for First-Order Conic Solvers) [13].We evaluate phase retrieval with and without COACS preprocessing based on TFOCS, realizing marked improvements.
In the remaining sections, we first describe our optimization problem in greater detail, followed by a detailed account of implementation concerns, experiments, results, and their evaluation.

Intensity healing using COACS
We will focus on the case of 2D imaging, 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 [14], and for that matter to further dimensions.
Assume that the 2D projection being imaged 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 center-cropped discretizations of their continuous 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.Figure 1 illustrates the effects of discretizing the infinite ideal diffraction pattern of our simulated test particle, with small but non-negligible artifacts permeating to the edge of the image.The physical processes that we model are based on eq. 1.To handle those properly in the discretized space, we extend eq. 2 with the Hann window, known from spectral techniques in signal processing.In 2D, the Hann window in real space amounts to a convolution with a 3x3 blurring stencil, and in the Fourier space a translated sinusoid curve dropping off to 0 at the edges.This ensures a representation where the high-frequency content smoothly goes to zero.Having no explicit window is equivalent to a sharp rectangular window, giving rise to the spectral leakage effects seen in Figure 1.
If we can assume a photon-counting detector with uniform quantum efficiency r, we will then have a sampled diffraction pattern in terms of an integer matrix B, where: where Po is a Poisson distribution with the rate parameter λ of appropriate dimension identifying the mean (and variance).For high rates, B can be treated as a Gaussian, or even as an exact observation B i,j = r|X i,j | 2 .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.

Support constraints
Since a single particle is imaged in isolation, we can introduce a support constraint, where P S , or P S , is 0 (where signifies the complement).For phase  retrieval, one can thus express this as: arg max Alternating-projection methods are very commonly used for phase retrieval.They were pioneered by Fienup [2], and are still popular in practical applications use, as e.g.implemented in [15].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, rather than stochastic as given above.Both of these constraints can be varied, to better explore the full solution space and accelerate convergence [2,3], or in order to account for measurement error in the intensities [16].More recent attempts include using the Alternating Direction Method of Multipliers for easily enforcing more general constraints in both spaces [17].
These attempts 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 as PhaseLift [10] and PhaseCut [12].These approaches claim to produce unique solutions under favorable conditions with a vanishing probability of not doing so, such as observing an object using alternating binary masks, or from unique angles randomly distributed over a sphere [10].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 pp, where p would be the flattened vector from our 2D matrix P. A megapixel image would thus give rise to a trillion-element matrix.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.A convex representation of that problem can be expressed as minimizing the trace.
The authors of PhaseLift note that while oversampling in theory gives a unique solution to the phasing problem, uniqueness is not equivalent to practical numerical attainability.In fact, they even claim, based on simulations including their own method as well as alternating projections [11] (original italics): The ill-posedness of the problem is evident from the disconnect between small residual error and large reconstruction error; that is to say, we fit the data very well and yet observe a large reconstruction error.Thus, in stark contrast to what is widely believed, our simulations indicate that oversampling by itself is not a viable strategy for phase retrieval even for nonnegative, real-valued images.
This argument is supported by additional claims about numerical stability.We will argue that a new analysis, of the support constraint specifically, will allow a much more favorable numerical treatment.

A support constraint for the Patterson function
The quantity that is sampled to form the directly observable B is Y = X X (the elementwise Hadamard product with the conjugate, forming the absolute value squared).In crystallographic applications, the Fourier transform of Y is referred to as the Patterson function.Using the convolution theorem, we can characterize this as: Here, * indicates convolution.The existing support constraint (4) for X can then be formulated as a (weaker) support constraint for Y: arg max Based on the structure of the probability distribution function in ( 3), ( 7) 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.The self-convolution effect of treating the intensity directly has also led some authors to refer to this as the autocorrelation function.

Implementation
If the noise structure would have been Gaussian, the resulting problem from (7) would have amounted to an ordinary least squares problem.A solution could 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.For a practical 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 would be impossible to transform back to amplitudes in the original phasing problem.
For a long time, so-called interior point methods have been a method of choice for many families of convex optimization problems [18], with free [19] as well as commercial [20] solvers, and more high-level modeling environments [21].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 would also have 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 instead chosen to rely on first-order methods, where only the gradient is computed.The TFOCS [13] package is a Matlab toolbox, and is presented by the author as a set of templates from which one can build tailored solvers.There are a couple of different canonical forms for formulating problems in TFOCS, and we have used the one that corresponds to the main function call tfocs: The linear operator A 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.h on the other hand, only has to be prox-capable, meaning one can (quickly) solve: In our context, f corresponds to the noise model, while h represents the support constraint.Both of these are implemented as custom functions, for numerical accuracy reasons outlined in the next section.

Numerical concerns
As noted by the authors in [11], 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 4 , where q is the scattering angle [22].The dynamic range needed for a problem in intensity space is higher than in the amplitude space, due to the quadratic relationship between the two.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 hold for the final result, but for properly evaluating the change undertaken between iterates.
The steps outlined below try to consistently reduce the chance of: • Truncation error due to very small terms being added to large terms.
• Loss of precision due to large negative values being added to positive values of almost equal magnitude, or vice versa.
• Very small steps due to the optimal solution being located at the border of the allowed domain.
The specific places where this has been taken into account can be summarized as: • Lipschitz backtracking structure in TFOCS • Translation of the solution variables and customized probability function • Relaxed quadratic barrier and accelerated continuation 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.

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 [13,23].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 backtracking is using a local, adaptive estimate of a Lipschitz constant for the gradient of the objective function.
As elaborated upon in [13], there are multiple ways to estimate a bound on the Lispchitz 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 that is 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.

Translation of the solution variables and customized probability function
The Fourier transform is a linear operator.Thus, the relationship between Y and P can be decomposed: In this case, the optimization variable for the first-order method in TFOCS can be Y .With Y 0 chosen properly, such as the resulting solution from a previous run of the TFOCS algorithm (the previous "outer iteration"), the norm of Y should be much smaller.This change reduces the effects of truncation errors in the gradient steps in individual iterations, including the probability of resorting to the conservative Lipschitz bound discussed above.It also reduces the effects of loss of precision due to large values of opposite sign canceling out within the Fourier transform.The effects do still exist on a global level, but they do not hamper the evaluation of objective values, constraint satisfaction, and search direction locally during the iteration process.The equivalent translation is not possible in the phasing problem, since any changes in the phase estimates affect the full amplitudes, whereas an additive decomposition is possible in the intensity space.
To fully realize the benefits of translation, the corresponding shift needs to be done in the computation of the f term of the objective function as well.We have implemented a custom log Poisson implementation, which computes the change in objective function value between Y 0 and Y using Y as input.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.
In our implementation, we have chosen to implement the windowing factor inside f , rather than in A. This makes the effects of the windowing on the resulting gradients more visible (the gradients are scaled by the window squared).
Since the reason to introduce l is to ensure a smooth behavior amenable to optimization, it makes sense that the effective gradient should be of equal magnitude everywhere.Therefore, the actual value used for the barrier width is inverted by the window, resulting in a more lenient barrier in the high-frequency areas, but also requiring an overall lower value of l to guarantee high accuracy at all angles.To stabilize the inverse process, the Hann window in amplitude space (which is squared in intensity space) has an extra term of 10 −3 in order to avoid singularities.

Relaxed quadratic barrier and accelerated continuation
The extra quadratic barrier 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).This value is then decreased by a factor of .5 repeatedly until reaching 2 −46 ≈ 10 −14 .Due to the Hann window scaling of l, this means an effective barrier of 10 −8 at the edge of the image.This scheme ensures a very quick convergence to the relaxed problem, followed by successive sharpenings 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 2−i ).Since the major trend influencing the results should be the reduced fluctuations around 0, the initial guess for iteration i, when i > 2, will be X i−1 + 0.5(X i−1 − X i−2 ), i.e. reproducing half the change that resulted from the previous halving of l.This acceleration scheme is similar to the ones more thoroughly investigated for the SCD model in [13].
Naturally, the actual X values used are translated 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 adaption.Therefore, a second set of inner continuations is applied after a fixed number of individual inner iterations of the optimization scheme.In these inner continuations, a separate acceleration scheme is used after the second iteration, where the new starting point is displaced with an accumulation factor of .9.If at any point the resulting X estimate is found to have a shorter distance to the old end point than the new starting point, the acceleration is deemed unsuccessful and that inner iteration is restarted without acceleration.We have found this ad-hoc scheme to give very reasonable results with an acceptable time usage.The inner continuation loop is terminated when the change in the objective function, scaled by the value of l, falls below a pre-defined threshold.

Window-aware support constraint
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.We therefore ventured into employing a strong, but non-infinite, quadratic penalty, chosen to be 5e7 l .
The quadratic penalty, as well as the original zero-projection constraint, has a drawback in terms of the spectral behavior.The sharp cutting at the edge of the support will induce high-frequency signals in the diffraction space, while we are already employing the Hann window to the original diffraction signal in order to reduce that signal.Therefore, we are using the same kind of twolevel Hann window intensity space for the support penalty matrix in the Fourier space, which amounts to a gradual decrease in support penalty over a width of 5 pixels.
In practice, similar results were achievable without this addition, but the time usage was much higher.This can be understood by the fact that in our simulations, the support was in fact wider than the object.Therefore, the ideal solution is already zero in the border region.For a fully converged, ideal, solution, the two constraint formulations are therefore identical.The windowed version will only accelerate convergence by avoiding spectral leakage.

Data and experiments
Simulated diffraction data were pre-treated 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 [24] package.The test particle chosen in the simulations consisted of an icosahedron, with a single very highdensity sphere placed on a point off-center on one of the vertices.The shortest diameter between vertices was 20 pixels, imaged on a virtual 256x256 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 the ability to resolve the feature at all.At the same time, our simulated sample is as a first approximation roughly spherical, 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 25x25 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 [25] 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 images.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.

Results
A qualitative comparison of the results possible using intensity healing are given in Figure 2, for the 50 patterns simulated.For each pair, the COACsed pattern is showed to the right.The high-density spherical feature on the edge is clearly visible in all healed reconstructions, and the contours of the icosahedron projection are reasonably traced, while the traditional phasing method is only able to produce an elongated shape without any detail, and a much more gradual fading into the outer areas of the support.
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 Figure 1.The phased COACS patterns and especially the nonphased COACS patterns are clearly superior to the patterns achieved using traditional phasing.There is no perceptible ringing effects from the highly regular rectangular support used in the COACS patterns.
The crystallographic R factor, or equivalently the L1 norm of the difference between the recovered wave amplitudes and the true amplitudes, normalized by the norm of the true amplitudes [26], is 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.In Figure 4, the R factors calculated at various radii (in pixels) are also shown.The traditional phasing method is 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 due to the shape of the high-density spherical feature in the pattern at 90 pixels, 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. The healing process required approximately 30 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.

Discussion
The time usage for the intensity healing is much longer compared to phasing, 30 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 5 ).It should be noted that the current COACS implementation is a straightforward one in Matlab using the TFOCS toolbox, while the phasing implementation libspimage [15] 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.The adaptive backtracking of the step-size in TFOCS does mean that the number of evaluations can be higher, but in our testing this was still lower than a total factor of 5. 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 CPUbased implementation in Matlab for intensity healing, and an efficient C++ implementation on GPUs using CUDA.Another avenue would be to use the existing port of TFOCS methodology into the TensorFlow framework.
We can also conclude that our current acceleration scheme is ad hoc rather than optimal.In fact, the translation of our problem into one of general convex optimization means that one can much more easily identify and benefit from existing acceleration schemes that have shown success in other applications.We have tried reformulations using the TFOCS SCD scheme with accompanying relaxations of µ (not shown), but those did not perform better than our ad hoc approach in terms of accuracy or performance.
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 Pois-  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.Peak at around 90 pixels is due to the R factor being a relative error measure.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.son probability distribution with a Gaussian component for values close to 0, 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 reconstruction processes that start to "walk" along the image.This is due to a non-ideal phase ramp being induced by the combination of the two constraints, resulting in movement in real space.This, in turn, means that the real space object will repeatedly "bump" into the support, with new artifacts being introduced as the real space constraint implementation tries to remove the signal.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 those.It is our opinion that a proper phasing method should be able to produce R factors similar to those we report for the COACS method alone in Table 1 and Figure 4.
The straightforward structure and generality of the convex optimization formalism as well as the TFOCS library in particular also make it feasible to add additional constraints.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.

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, which 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 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.

Figure 1 :
Figure 1: Simulated particle with and without Hann windowing.a) Original simulated diffraction pattern with high-density sphere.b) Center crop 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 windowing applied, resulting in lower intensity away from the center.e) Discrete Fourier transform of pattern after windowing.Slight blur visible.f) The image in e) with limited range.Artifacts found in c) mostly absent.

Figure 2 :Figure 3 :
Figure 2: 50 phased reconstructions of sparse patterns based on the simulated particle.Pairs of phased healed pattern using COACS (left) and phased original patterns (right).Average of 10 best individual phasings out of 100 replicates for each of these.

Figure 4 :
Figure4: 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.Peak at around 90 pixels is due to the R factor being a relative error measure.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.

Table 1 :
R factor means and standard deviations.