Facing the phase problem in Coherent Diffractive Imaging via Memetic Algorithms
Alessandro Colombo Davide Emilio Galli Liberato De Caro Francesco Scattarella Elvio Carlino




alessandro.colombo6@unimi.it
Università degli Studi di Milano, via Giovanni Celoria 16, 20133 Milano, Italy


davide.galli@unimi.it 
Università degli Studi di Milano, via Giovanni Celoria 16, 20133 Milano, Italy


liberato.decaro@ic.cnr.it
Istituto di Cristallografia, Consiglio Nazionale delle Ricerche  (IC-CNR), via Giovanni Amendola, 122/O, 70126 Bari, Italy


scattarella@iom.cnr.it
Istituto Officina dei Materiali, Consiglio Nazionale delle Ricerche (IOM-TASC-CNR), Strada Statale 14 km 163.5, 34149 Trieste, Italy 


carlino@iom.cnr.it
Istituto Officina dei Materiali, Consiglio Nazionale delle Ricerche (IOM-TASC-CNR), Strada Statale 14 km 163.5, 34149 Trieste, Italy




Coherent Diffractive Imaging is a lensless technique that allows imaging of matter at a spatial resolution not limited by lens aberrations. 
This technique exploits the measured diffraction pattern of a coherent beam scattered by periodic and non–periodic objects to retrieve spatial 
information. The diffracted intensity, for weak–scattering objects, is proportional to the modulus of the Fourier Transform of the object scattering function. Any phase information, needed to retrieve its scattering function, has to be retrieved by means of suitable algorithms. 
Here we present a new approach, based on a memetic algorithm, i.e. a hybrid genetic algorithm, to face the phase problem, which exploits 
the synergy of deterministic and stochastic optimization methods. The new approach has been tested on simulated data and applied 
to the phasing of transmission electron microscopy coherent electron diffraction data of a SrTiO_3 sample. 
We have been able to quantitatively retrieve the projected atomic potential, and also image the oxygen columns, which are not directly visible 
in the relevant high-resolution transmission electron microscopy images. Our approach proves to be a new powerful tool for the study of matter 
at atomic resolution and opens new perspectives in those applications in which effective phase retrieval is necessary.


Facing the phase problem in Coherent Diffractive Imaging via Memetic Algorithms
    Elvio Carlino
    December 30, 2023
===============================================================================









§ INTRODUCTION

Full-field and scanning microscopes can be either lens-based or  lensless imaging systems. Coherent Diffractive Imaging (CDI) is a  lensless technique that  permits imaging matter at a spatial resolution not limited by lens aberrations. The seminal idea of CDI was due to David Sayre in 1952 <cit.> but it was only experimentally demonstrated for X-rays in 1999 <cit.> and, more recently, also for electrons, using a Transmission Electron Microscope (TEM), giving rise to the Electron Diffractive Imaging (EDI) <cit.>. 

The goal is to retrieve a qualitative/quantitative image of a scattering function related to a physical property of the scattering object, such as the electron density (X-ray CDI) or the atomic potential (EDI).
High Resolution TEM (HRTEM) images of the projected atomic potential are phase-contrast images limited by the high-order aberrations of the objective lens, which distort the phase of the scattered wave function, giving rise to images of the sample, which in general are not immediately interpretable in terms of its atomic structure <cit.>.
Instead, diffraction patterns of scattering objects are not affected by these aberrations. Therefore they contain, in principle, undistorted information on the scattering function at a better spatial resolution with respect to lens-based imaging systems <cit.>. The diffracted intensity, for weak-scattering objects, is proportional only to the modulus of the Fourier Transform (FT) of the scattering function. Any phase information, which is experimentally lost (phase problem <cit.>), has to be retrieved by means of suitable algorithms. The lensless image of the scattering function, obtained by means of an inverse FT of the diffraction pattern once that the correct phase has been retrieved, is characterized by a final resolution experimentally limited only by the Numerical Aperture (NA_diff) corresponding to the highest spatial frequency contained in the diffraction pattern that can be related to the atomic structure of the investigated sample <cit.>. Wavelength, noise, radiation damage, thermal and mechanical stability of the experimental setup, dynamics of the detector, etc. <cit.> could limit the spatial resolution achievable. 

In order to find a unique solution to the phase problem, phase retrieval algorithms need a-priori constraints, such as fixing a region around the sample characterized by zero scattering <cit.>. The extension of the zero-scattering region has to be large enough regarding the object support S, containing the scattering object. This requirement is the so-called oversampling of the diffraction pattern, needed to satisfy the Nyquist sampling requirements <cit.>.


For extended samples, alternative approaches have been developed to have enough a-priori information to make the phasing problem overdetermined. In particular, Abbey et al.   <cit.> for X-rays CDI experiments substituted the region of zero scattering with a region of zero illumination using a confined probe with a technique named keyhole CDI. The possibility to perform Keyhole coherent diffraction experiments has been demonstrated also for electrons in EDI, obtaining the KEDI approach <cit.>.

So far, several phase retrieval algorithms have been developed <cit.>, which are mostly evolutions of Fienup’s modification <cit.> of the Gerchberg-Saxton's algorithm <cit.> working with a dual-space strategy. 
The common approach of many phasing algorithms is to impose constraints both in real space (the prior information on the zero scattering/illuminating region) and in Fourier space (the amplitudes are adapted to the experimental values). The imposition of a constraint in one space always causes the violation of the constraint in the other. Consequently, the standard strategy is to use iterative schemes but, in this way, the global minimum of the reconstruction errors is often reached with difficulty <cit.>. Indeed, the support S of the scattering function is either unknown (X-ray CDI) or known at a worse spatial resolution. This is the case of EDI, for example, where the support is obtained by means of a lens-based image of the scattering function (projected atomic potential), i.e. by the HRTEM image. 

In case of unknown complex scattering functions the knowledge of the support S (non-zero illumination region) at the same spatial resolution of the measured diffraction pattern seems to be mandatory for the success of the phasing, especially when phasing algorithms are not suitably structured to escape stagnation in local minima <cit.>. This aspect would limit the advantages of indirect imaging based on lensless systems with respect to conventional lens-based set-ups.

Indeed, standard phasing approaches, such as the Hybrid Input-Output (HIO), are mainly deterministic iterative algorithms, which, going back and forth from the real to the Fourier space, try to optimize a specific error functional <cit.>. For this reason they are highly efficient in finding local minima, but they suffer  from stagnation mainly due to the incomplete knowledge of the support S. Furthermore, the final result is highly dependent on the initial conditions <cit.>.
In order to overcome these limitations, phasing procedures usually involve a lot of parallel and independent retrieval processes, with different initial conditions, choosing the scattering function with the lowest reconstruction errors as a possible solution.
A first step toward a smarter use of information coming from multiple phase retrieval processes is the Guided Hybrid Input-Output algorithm <cit.>.

A different approach to face the phase problem in CDI could make use of pure stochastic optimization methods. However, a major limitation of these methods is their inefficency once the number of unknowns is huge, so that, in typical CDI applications, such an approach is doomed to fail even by exploiting actual super–computing facilities.

In this article we make a step further proposing a new hybrid stochastic approach to better explore the phase solution space through a smart use of Genetic Algorithms (GAs) <cit.>. GAs have been already applied to the phase problem in different fields <cit.>.
The novelty of our new approach consists in the development of a Memetic Algorithm (MA) <cit.> 
in the context of phase retrieval applied to CDI; this scheme represents a natural choice for a smart merging 
of stochastic and deterministic optimization methods: the algorithm has been developed hybridizing a GA, which guarantees a wide exploration of the configuration space, with local optimization algorithms like Hybrid Input-Output and Error Reduction.

We have shown on simulated data that the MA phasing approach is able to retrieve the correct scattering function, when it is real, imposing a very loose support constraint S in the direct space. 
Moreover, still on simulated data, we have shown that the new MA phasing approach is able to retrieve the correct scattering function, even if it is complex, starting from a knowledge of the support S at a resolution four times worse than the one corresponding to NA_diff, even worse than normally observed in EDI/KEDI real experiments. Our new approach shows convergence performances towards the global minimum that go well beyond those achievable by standard phasing deterministic algorithms. 

Finally, we have applied the GA-based phasing approach to a KEDI experiment realized on a SrTiO_3 sample in a [100] axis orientation.
The image obtained after the phase reconstruction is a detailed structural map of the specimen atomic potential projected along the [100] direction  at a sub-Ångström spatial resolution corresponding to the highest frequency measured in the experimental diffraction pattern. The intensity distribution enables one to distinguish between atomic sites containing different chemical species. Also the oxygen signals can be detected, despite the presence of heavy atoms in the crystal cell, whereas they are not visible in the relevant HRTEM image.

These results pave the way to the highest spatial resolution, accuracy and reliability achievable in lensless imaging and represent a new powerful tool for the study of matter.



§ RESULTS AND DISCUSSION




 §.§ The Memetic Phase Retrieval approach 


A GA <cit.> is a stochastic optimization method that imitates the survival–to–fitness typical of the natural evolution of a population. In general, this is obtained by elaborating the genetic information via three genetic operators: Selection, Crossover and Mutation.
In our algorithm we induce the genetic dynamics on a set of initial densities {ρ_i(x⃗) }_i=1 … N_p (also called population) which represents N_p possible solutions to the phase problem. 
In standard GAs Mutation and Crossover induce a stochastic shift on every element of the population in the space of configurations; this improves the ability of exploring the space, but makes the GA efficient only when N_p increases with the number of unknowns <cit.>.
In the phase retrieval problem this condition makes standard GA impracticable due to the computational cost.

A practised way to overcome this issue consists in implementing hybrid GAs, known as Memetic Algorithms (MAs) <cit.>, where global optimization is boosted by local optimization procedures. This operation, in the framework of Memetic Algorithms, is known as Self (or Local) Improvement.



Using this method, we introduce local iterative deterministic phase retrieval procedures in our algorithm as Self Improvement operation.
For this reason, the new proposed phasing approach is called Memetic Phase Retrieval (MPR). Fig.<ref> shows an overview on the procedure, while a more detailed description of the method is reported in the Methods section. It is worth noting that the standard phasing approach can be seen as MPR without the genetic operators Selection, Crossover and Mutation. 
MPR has to be considered as a smart framework to take advantage and improve the performances of any current iterative phase retrieval approach; it is clear, in fact, that any phase retrieval algorithm can be implemented in MPR.
In this work we use the Error Reduction and Hybrid Input-Output algorithms as Self Improvement operations because they are simple, well known and well-characterized methods. Moreover, as we will show in the following, the sole inclusion of HIO and ER inside MPR is enough to build a very powerful phase retrieval algorithm. 
MPR actually includes also methods for the retrieval of the optimal support function, like the Shrinkwrap algorithm <cit.>; this feature makes MPR a Co-evolving Memetic Algorithm <cit.>, where the Self Improvement co-evolves along with candidate solutions. This peculiar feature will be discussed in future works because it has not been used in the present application on KEDI data, where sufficient information about the support function was available, making Shrinkwrap procedure unnecessary.







 §.§ Testing MPR on simulated data

In the Supplementary Information we discuss in detail several numerical tests performed to verify the potentialities of the new proposed approach to reliably retrieve in a reliable way phase information in comparison to standard deterministic phasing approaches. Here, we reassume the main results obtained by simulations before discussing the application of MPR on true experimental data of a KEDI experiment.

In the comparison of the performances between MPR and standard phasing algorithms, applied on simulated data, we will focus on the best phase reconstructed by both methods and not on a statistical analysis of the results obtained from the set of phase retrieval processes started in parallel with different initial conditions. This is necessary because the genetic dynamics has the distinctive feature to mix and share parallel information during the stochastic evolution; this makes smarter the stochastic search for a better phase, but also tends to push the population near to the best candidate solution.
A statistical analysis would be thus biased in favor of MPR. 
Moreover, we are going to compare MPR and standard phasing algorithms under equal conditions: the same number of candidate solutions in the population and the same kind of iterative phase retrieval algorithms.





The first test of comparison between MPR and the standard phasing approach is a phase retrieval of a positive and real-valued two-dimensional (2D) scattering function, as shown in Fig.<ref>a.
In the test the provided support function S (Fig.<ref>b) is a square four times smaller than the total area of the direct space, which gives a constraint ratio <cit.> Ω=2.
With this regard it is useful to remember that Ω≥ 1 is the mathematical condition to assure the existence of a unique solution to the phase problem. Ω is defined as the ratio between the support of the autocorrelation of the scattering object and two times the object support, which are areas for 2D phase problem, volumes for 3D ones. In this example, the support function is not updated during the phase retrieval.





Fig.<ref> shows the obtained results for a population N_p = 16384. In particular, Fig.<ref>b shows the support S, never updated during the phasing process. Fig.<ref>a shows the retrieved unknown function obtained by standard phasing algorithms (see Supplementary Information for further details). The poor quality of the retrieved phase is due to the weak constraint in real space, as the support constraint S has not been updated. Fig.<ref>c shows the retrieved unknown function obtained by MPR (see Phase retrieval of real-valued data in Supplementary Information for further details). 


The performances of the new proposed stochastic approach are much better than the classic deterministic phasing methods. MPR works accurately even without a tight support constraint.

In the case of simulated data one can evaluate the true error defined as the normalized absolute difference for every pixel of the retrieved 2D function with respect to the solution (Lena image, adapted from the picture 4.2.04 in the USC–SIPI image database<cit.>) (see Supplementary information for more details).
The true error for the deterministic phasing approach is larger than 24%. Instead, MPR leads to a true error less than 1%.






A second numerical test of MPR concerns the reconstruction of a complex-valued scattering function.
In particular, we have considered a situation typically encountered in EDI, in which HRTEM allows us to obtain a lens-based image of the sample under investigation characterized by a worse spatial resolution with respect to that corresponding to the NA_diff of the measured diffraction pattern. In order to simulate this experimental situation, we have binned with a factor 4 the module of the scattering function to be retrieved (Fig.<ref>a) to obtain a rough estimation of its support, thresholding the binned image as shown in Fig.<ref>c, whereas Fig.<ref>b shows the phase in direct space that has also to be retrieved. 



This data can be represented as depicted in Fig.<ref>a via the Hue-Saturation-Value (HSV) color system, where the information on the phase is stored in the hue, the modulus corresponds to the value and the saturation level is set to the maximum.
Even in this case, the standard phasing approach is far from recovering the correct complex scattering function, reported in Fig.<ref>b. Instead, MPR is able to correctly retrieve both module and phase of the complex unknown scattering function (Fig.<ref>c) (see Supplementary Information for further details).
The final true error for this test is about 10% for the standard approach and 1.5% for MPR.







 §.§ Application of MPR on experimental data for Keyhole Electron Diffractive Imaging


KEDI experiments <cit.> are challenging tests for phase retrieval algorithms as the scattering function to be reconstructed is complex. SrTiO_3 was considered as a case study for the great importance of this oxide from both an applicative and a fundamental point of view. The role of the oxygen sub-lattice is of particular importance in the studies of two dimensional electron gases formed at the interface between two insulating oxides and has recently attracted great attention <cit.>. The capability to image the lattice of complex oxides at atomic resolution is necessary to understand the intriguing properties of this class of material. Moreover imaging of a light chemical element, such as oxygen, in a matrix of heavier atoms, like titanium and strontium, is not straightforward <cit.>. The samples were prepared for KEDI experiments in a [100] zone axis as this configuration enables the imaging of different atomic species in the crystal sub-lattice (see Methods section).







KEDI requires an HRTEM image (Fig.<ref>a) and a nano-diffraction pattern acquired from the same sample area with the same electron optical conditions <cit.>.
The HRTEM experiment enables us to image the sample (Fig.<ref>a and <ref>b), to complement the diffraction pattern at the lower spatial frequency (Fig.<ref>c) and to estimate the support (Fig.<ref>d) at the resolution allowed by the experimental conditions and by the electron objective lens aberrations <cit.>. In the case of the electron-optical set-up used for these experiments the relevant spatial resolution in the HRTEM image at optimum defocus is 0.19 nm <cit.>. The HRTEM image in Fig.<ref>b has been successfully simulated in the framework of full dynamical Bloch-wave approach <cit.> for a thickness of 25 nm and an underfocus value of 41.3 nm (see Supplementary Information). It should be noted that the phase contrast in the HRTEM image of Fig.<ref>b does not show any evident clues that could be correlated to the presence of the oxygen atomic columns which should be seen in the [100] projection of the SrTiO_3 atomic potential, as evidenced in the simulation (see Fig. S7 of Supplementary Information).
Indeed, it is worthwhile to remark that the HRTEM image is, in general, an interference pattern of the waves scattered by the atomic potentials in the specimen, and the positions of the maxima and minima in the image cannot be straightforwardly interpreted as structural features and therefore the comparison with the simulated images is needed <cit.> (see Fig. S7 of Supplementary Information).

The KEDI diffraction pattern, shown in Fig.<ref>c has been obtained by combining the measured diffraction pattern with the modulus of the HRTEM image FT, after a suitable matching procedure requiring its rotation and scaling <cit.>. 
The pattern in Fig.<ref>c is the starting point for the phase retrieval process. 
It is worth noting that the MPR phasing process has been carried out without any a–priori information about the phases, an information which is, instead, needed by standard phasing procedures applied to KEDI <cit.>.
Further details on MPR applied to experimental data have been reported in the Supplementary Information.



Fig.<ref> shows the retrieved scattering function obtained by using MPR, where the brightness corresponds to the modulus and the hue to the phase of the retrieved real-space complex-valued scattering function. 
The long range phase variation is due to the phase variation of the illumination nano-probe <cit.>.




Fig.<ref>a shows the phase retrieved amplitude for the structure of the SrTiO_3 seen in a [100] projection.
Fig.<ref>d has been obtained by subtracting the contribution of the TEM illumination function to quantify the SrTiO_3 projected potential.
The first important point that should be emphasised, is that in the phase recovered image the positions of the maxima are correctly in correspondence with the expected positions of the atomic columns seen in the [100] projection. In other words, the phase reconstructed image is a structural image of the specimen. This is a fundamental issue that paves the way for quantitative structural imaging at atomic resolution. 
Indeed, as shown in Fig.<ref>d, the Sr and Ti+O columns are precisely seated on the relevant square sublattice of the SrTiO_3 in the [100] projection (see Fig.<ref>c and Fig.<ref>b). Approximately in the center of the sublattice there is a lower signal which corresponds to the oxygen columns.
A second point concerns the retrieval of quantitative information about the atomic potential. By comparing data shown in Fig.<ref>d with the expected projected atomic potential (Fig.<ref>c) we found that the ratios between the intensities of Sr, Ti+O and O columns are correctly retrieved, providing truly quantitative information on the specimen. In particular, the expected intensity ratios are I_O / I_Sr=0.35 and I_Ti+O / I_Sr=0.96  while the experimental retrieved data give I_O / I_Sr=0.35±0.05 and I_Ti+O / I_Sr=0.89±0.10.
The possibility to do KEDI experiments by MPR reconstruction paves the way for a detailed structural characterization of the investigated samples and opens up new possibilities for the understanding of the properties of the matter at sub-Ångström resolution.






 §.§ Conclusions

In this work we have discussed some of the potentialities of a new phasing approach, stochastic in the exploration of the space of solutions, based on a memetic algorithm, applicable both to X-ray and electron coherent diffraction imaging. We have tested the new phasing algorithm, named Memetic Phase Retrieval, on simulated data and we have obtained the result that the knowledge of the support - which defines the boundaries in the direct space of the unknown scattering function that one wants to retrieve - is less binding than previously reported. The more efficient exploration of the space of solution, possible thanks to the stochastic genetic procedures implemented in MPR, should be the higher gear of the new proposed phasing method with respect to those already available. Indeed, both the possibility to correctly retrieve a real-valued scattering function by its diffraction pattern without imposing any tight support and to reconstruct a complex-valued scattering function by using a low-resolution estimate of the support, are examples of the great capabilities of MPR to face the phase problem. Our tests on simulated data demonstrate the superior capabilities of the MPR for accurate phase retrievals. Indeed, by using the same computational resources, the MPR approach has proved to be much more powerful than deterministic phasing procedures in facing the phase problem. 
The application to an experimental case of Keyhole Electron Diffraction Imaging has shown that the atomic potentials of SrTiO_3 can be quantitatively imaged, representing a relevant improvement for the study of the matter.
We believe that the Memetic Phase Retrieval approach could be of interest in all the fields that require accurate phase retrievals.




§ METHODS




 §.§ The phase problem as an optimization problem 

The ideal solution to the phase problem is a function ρ_s (x⃗), representing the spatial distribution of the sample, whose Fourier Transform (FT), ρ̃_s (q⃗), has a square modulus equal to I(q⃗), which is proportional to the experimental diffraction pattern intensity. 
ρ_s (x⃗) is also assumed to be zero outside a well-defined region of the real space, the so-called support S, in order to satisfy the oversampling condition, which assures the necessary information to retrieve ρ_s (x⃗) <cit.>. 
Here, x⃗=(x,y), with x and y the cartesian components of the position vector x with respect to the reference system. 
Analogously, q⃗=(u,v), where u and v are the spatial frequencies components with respect to the reference axes.

In principle, this solution can be represented as an intersection of sets.
ℳ can be defined as the set of all functions ρ (x⃗) compatible with experimental data I(q⃗), i.e:

    ℳ = {ρ (x⃗) : | ρ̃ (q⃗) |^2 = I(q⃗) }.



𝒮 is, instead, the set of all functions ρ_s (x⃗) satisfying the oversampling condition, which is defined by a binary function Π (x⃗) representing the object support S. So, the set 𝒮 is described by:


    𝒮 = {ρ (x⃗) : ρ(x⃗)= Π (x⃗) ρ(x⃗) }.


Thanks to (<ref>) and (<ref>), the ideal solution is:

    ρ_s (x⃗) = ℳ∩𝒮.


The main issue concerning experimental measurements is the presence of noise and lack of data; this, in general, implies that 𝒮 and ℳ do not intersect:
 
    ℳ∩𝒮 = ∅.



Due to the condition (<ref>) a different way to define what we mean by “solution” is needed. 
It is useful, at this point, to introduce two projection operators, which act on the function ρ (x⃗):


    P_ℳ : P_ℳρ(x⃗) = ℱ^-1 [ √(I(q⃗)) e^i [ρ̃ (q⃗)]](x⃗),



    P_𝒮 : P_𝒮ρ(x⃗) = Π(x⃗) ρ(x⃗).


It's trivial to prove that P_ℳ and P_𝒮 are projectors on sets ℳ and 𝒮, previously defined in (<ref>) and (<ref>).
Thanks to these operators, it's now possible to give a new definition of solution in place of the one defined in the eq. (<ref>):

    ρ_s (x⃗) = min_ρ D[P_ℳρ(x⃗), P_𝒮ρ(x⃗)],

where the functional D[A,B] represents the metric of the space. Hereafter, we will refer to the eq. (<ref>) whenever we will talk about the “solution” of the problem, ρ_s (x⃗).

It's now clear that, in this framework, finding a solution to the phase problem means minimizing the distance between sets ℳ and 𝒮: the phase problem becomes an optimization problem for the quantity D[P_ℳρ(x⃗), P_𝒮ρ(x⃗)], which can be reinterpreted as the error of the recovered density ρ(x⃗).
Different definitions of the metric imply different definitions of the error assigned to a given ρ(x⃗) and, as consequence, different optimization targets.
We can define the error functional E[ρ] as

    E[ρ] = D[P_ℳρ(x⃗), P_𝒮ρ(x⃗)],

such that the eq. (<ref>) turns into

    ρ_s (x⃗) = min_ρ E[ρ].


Standard approaches to the phase problem are mainly deterministic iterative algorithms which, going back and forth from the real to the Fourier space, try to minimize a specific error functional <cit.>. These methods are highly efficient in finding local minima, but they suffer from stagnation and the final result is highly dependent on the initial conditions <cit.>.
In order to overcome these issues, phasing procedures usually involve a lot of parallel and independent retrieval processes with different initial conditions and then selecting the one with the lowest error.

The founding idea of the new proposed phasing method is to better perform this parallel exploration of the space, through the use of a Memetic Algorithm.





 §.§ Selection as a Rigged Roulette

The Selection process is a delicate step in the Evolution process. A Selection strongly favoring only the better elements in {ρ_i(x⃗) }_i=1 … N_p (i.e., elements with the better fitness value) will improve the convergence speed, but the algorithm will suffer with stagnation in local minima.
On the other side, a selection process that weakly favors those elements will have, instead, an unstable convergence and will require an excessive length of time to find the solution.

There are several ways to select elements depending on their fitness value. The one chosen in this work is the so-called “rigged roulette”.
Once an error value E_i is assigned to every ρ_i(x⃗) in {ρ_i(x⃗) }_i=1 … N_p according to the eq.(<ref>), the set {ρ_i(x⃗) }_i=1 … N_p is ordered by increasing values of E_i (which is equivalent to a decreasing values of the fitness).
Whenever the algorithm has to select a ρ_i(x⃗) in {ρ_i(x⃗) }_i=1 … N_p for the Crossover operation, an index is extracted through the relation

    s = ⌊{rand[0,1)}^r · N_p ⌋ + 1 ,

where r ≥ 1 is related to the “strength” of the selection process. Usual values of r range from 1.5 to 2.5.

Eq. (<ref>) maps a flat distribution in [1, 0) ⊂ℝ to an unbalanced distribution in { 1, N_p }⊂ℤ, where the higher the value is the greater the probability is of getting a lower index and, therefore, selecting a better element.




 §.§ Differential Crossover

In the Natural Evolution process, the Crossover operation is the mixing of the parents' genetic pool.
In our implementation chromosomes are represented by every single (complex) value of ρ̃ (q⃗) = ℱ[ρ] (q⃗).
This means that, given two parent functions ρ_1(x⃗) and ρ_2(x⃗) selected according to their fitness, the son function ρ_son(x⃗) is created according to


    ρ̃_son(q⃗) = 
    ρ̃_1 (q⃗),    [0,1) > C 
    ρ̃_2 (q⃗),    ,

where rand[0,1) is a random number with flat distribution in [0,1) and C is a balancing coefficient between 0 and 1.

An improvement in performances can be obtained using the so called Differential Crossover <cit.> where, instead of selecting two parents, four parents, ρ_1(x⃗) ρ_2(x⃗) ρ_3(x⃗) and ρ_4(x⃗), are chosen.
The differential crossover acts as follows:


    ρ̃_son(q⃗) = 
    ρ̃_1 (q⃗),    [0,1) > C 
    ρ̃_2 (q⃗) + 
       D_c · [ρ̃_3 (q⃗)-ρ̃_4 (q⃗)]    ,

where D_c is called differential coefficient with typical values between 0.5 and 1.5.

The population of sons can be, in general, smaller than the whole population. This means that if the population of sons has N_s = G · N_p elements, a fraction of N_p - N_s parents, chosen randomly, will survive to the next generation.

The parameter G, which can be called genetic fraction, has values between 0 and 1: G=1 means that all of the parent population is replaced by the sons, while G=0 means that no sons are created, the genetic operators are switched off and we get a situation equivalent to the standard deterministic approach, as depicted in Fig.<ref>. 



 §.§ Mutation


Every element in the population {ρ_i(x⃗) }_i=1 … N_p may be subjected to a stochastic modification.
In this work, the mutation operation has been switched off because it does not introduce a remarkable improvement in the performance of MPR on treated data. Different implementations of the mutation operator are under study and will be topics of future works. 





 §.§ Self improvement via deterministic optimization
 
Optimization algorithms such as Error Reduction (ER) and Hybrid Input-Output (HIO) are efficient methods to find local minima or, more precisely, minima bounded to a region of the configuration space near the starting point. 
These algorithms are strictly bounded to the metric D[A(x⃗), B(x⃗)] defined as:

    D[A(x⃗), B(x⃗)] = ∫ dq⃗  [  |ℱ[A](q⃗)| -  |ℱ[B](q⃗)| ]^2.


This implies that the local optimization target is the functional E[ρ] defined as:

    E[ρ] = D[P_ℳρ(x⃗), P_𝒮ρ(x⃗)] = 
     = ∫ dq⃗  [ √(I(q⃗))  -  |ℱ[Πρ](q⃗)| ]^2.


In our algorithm this local optimization is carried on by elaborating every ρ_i(x⃗) with N_HIO iterations of the Hybrid Input-Output algorithm and N_ER iterations of the Error Reduction algorithm. 
In this work the global optimization target, i.e., the fitness of MPR, coincides with the local optimization target of ER and HIO algorithms just shown in (<ref>). 
This is not to be taken for granted because, in general, we can define any arbitrary global optimization target different from the local one (<ref>).
We are testing different fitness definitions for the global optimization, like the Csiszar's Information Divergence<cit.>,  and different local optimization algorithms.




 §.§ The choice of the initial guess

The phase retrieval process can be divided into two main steps. 
The first one concerns the choice of the initial population of densities {ρ_i(x⃗) }_i=1 … N_p. In the second step, we have to choose the parameters both of the genetic and the local optimization algorithms.

Standard approaches like Hybrid Input Output and Error Reduction need a single initial guess, which represent the first estimation of the solution.

MPR approach requires, instead, a set of initial guesses.
This set is produced from a single guess, simply randomly shifting every phase.
This means that, given an initial guess ρ_init (x⃗), every  ρ_i(x⃗) in {ρ_i(x⃗) }_i=1 … N_p is created via the relation

    ρ̃_i(q⃗_j) = √(I(q⃗_j))exp(√(-1)ϕ_j)  with ϕ_j = 
     =[ρ̃_init(q⃗_j)] + R_c ·rand[-π, π].


The parameter R_c, which has values between 0 and 1, depends on the accuracy of the initial guess ρ_init (x⃗).
If ρ_init (x⃗) is already a good estimation of the solution, it will be useful to set a low value (usually near to 0) for coefficient R_c in order to well explore the space near ρ_init (x⃗).
If, instead, ρ_init (x⃗) is considered to be far from the solution, it is useful to set a value of R_c near to 1, in order to explore also areas of the space far from ρ_init (x⃗).






 §.§ KEDI experiment


A KEDI experiment was performed following the procedure reported in <cit.> , which enables one to deliver a low dose of electrons to the specimen. The experiment requires the acquisition of an HRTEM image and a diffraction pattern by using the same electron optical set up. The experiments were performed by using a JEOL JEM 2010 F UHR operated at 200 kV. The cathode is a high coherence Shottky type. The microscope has an objective lens with low spherical aberration coefficient C_s=0.47 ± 0.01 mm and a relevant resolution at optimum defocus in HRTEM of 190 pm. The environment around the microscope is thermally and mechanically very stable allowing us to achieve in the scanning TEM (STEM) high angle annular dark field (HAADF) mode a resolution of 126 pm, which is the theoretical limit for the used electron optical set up <cit.>. 

In a KEDI experiment the optical setup produces an electron nano-beam. The latter defines the mathematical support of the scattering function for the illuminated nanometric region of the extended crystal. As in a microscope the field of view is proportional to the inverse of magnification, the size of the illumination function (beam size) is somehow related to the spatial resolution. The electron beam size S (which defines the support) is directly related to the final resolution to be achieved and to the size of the detector used to record HRTEM image and n-ED pattern. In fact, if the highest frequency of the diffraction signal recorded in the reciprocal space is ρ^-1pm^-1, we should have ρ at least two or three times the pixel size Δ_map of the phased map to have an electron projected potential two-dimensional map calculated with a sufficient number of points to be plotted continuously. 

For example, if we reached a final resolution – after the phase retrieval process – of ρ=70 pm we should have Δ_map∼ 25 ÷ 30 pm which, multiplied by the detector pixel number along a line, N=1024 in our experimental case, would lead to a spatial region O (scattering region plus non-illuminated surrounding region) of ∼ 25 ÷ 30 nm in size. Moreover, for the Nyquist theorem’s requirement, the illuminated beam size S (the support) has to be less than 2^-1/2 O, i.e., at maximum ∼ 17 ÷ 20 nm in size. Hence, in order to properly run the phase retrieval algorithms, the illuminated region of the sample in the direct HRTEM image has to be properly chosen with respect to the whole detector area to satisfy the above KEDI oversampling condition. 
Here, the cathode emission condition and the electron optical illumination system of the microscope has been experimentally set up to increase the probe coherence on the smallest illuminated area achievable <cit.>. The microscope has an illumination system composed by three magneto-static lenses. These lenses were operated independently, together with the electrostatic lens of the emitter, to produce the smallest-sized probe on the focal plane of the pre-field of the objective lens and hence the smallest-sized coherent parallel beam on the specimen. The emission conditions of the microscope cathode were chosen to increase the coherence of the electron probe by decreasing the temperature of the emitting tip. We used a heating current for the filament that halves the emission current with respect to the standard operation, decreasing at the same time the electron dose delivered to the specimen. The current density on the specimen was below the detection limit of the amperometer connected to the phosphorus screen of the microscope ( < 0.1 pA cm^-2), allowing us to acquire the relevant diffraction pattern on the 1024x1024 Charge-Coupled Device (CCD) camera without using the beam stopper for the direct beam. Thus all the diffracted intensities were available for the phasing process and a very small dose is delivered to the specimen. The small electron probe, without any changes, was used to acquire both HRTEM image and diffraction from the same area of the specimen. Fig.<ref>a shows the HRTEM image. The illuminated area is 10±2nm. The interference pattern of the phase contrast HRTEM image formed in the image plane of the objective lens is shown at a higher magnification in Fig.<ref>b In Fig.<ref>c the diffraction pattern formed in the back focal plane of the objective lens is shown. The central part of the pattern has been replaced, after proper scaling and rotation, by the FFT of the HRTEM image in Fig.<ref>a, as established in EDI method  <cit.>. The highest Miller’s index spot measurable in the pattern is the (5,5,0), which corresponds to a spacing of 55 pm. Thus, the expected gain in resolution of the maximum spatial frequency contained in the diffraction pattern (∼NA_diff^-1) with respect to that corresponding to the FT of the HRTEM image is about four times.






§ ACKNOWLEDGMENTS

This work was supported by the NOXSS PRIN (2012Z3N9R9) project and Progetto premiale MIUR 2013 USCEF.
We acknowledge the CINECA and Regione Lombardia LISA award LI05p-PUMAS, 
for the availability of high-performance computing resources and support.



§ AUTHOR CONTRIBUTIONS STATEMENT


A.C. and D.E.G. conceived the algorithm. A.C., supervised by D.E.G. and L.D., developed, implemented and tested the algorithm, and carried out the reconstructions. E.C. designed and performed the TEM EDI experiments, contributing to the optimization of the relevant data reduction and phasing. L.D. and F.S. made data reduction. All authors have equally contributed to the preparation and the revision of the text.



§ ADDITIONAL INFORMATION

Competing financial interests: The authors declare no competing financial interests.

How to cite this article: Colombo, A. et al. Facing the phase problem in Coherent Diffractive Imaging via Memetic Algorithms. Sci. Rep. 7, 42236; doi:10.1038/srep42236 (2017).







10
url<#>1urlprefixURL 


Sayre1952
authorSayre, D.
titleSome implications of a theorem due to Shannon.
journalActa Crystallographica
  volume5, pages843 (year1952).
<http://dx.doi.org/10.1107/S0365110X52002276>.

Miao1999
authorMiao, J., authorCharalambous, P.,
  authorKirz, J. & authorSayre, D.
titleExtending the methodology of x-ray crystallography to
  allow imaging of micrometre-sized non-crystalline specimens.
journalNature volume400,
  pages342–344 (year1999).
<http://dx.doi.org/10.1038/22498>.

Zuo2003
authorZuo, J., authorVartanyants, I.,
  authorGao, M., authorZhang, R. &
  authorNagahara, L.
titleAtomic resolution imaging of a carbon nanotube from
  diffraction intensities.
journalScience volume300,
  pages1419–1421 (year2003).

Huang2009
authorHuang, W., authorZuo, J., authorJiang,
  B., authorKwon, K. & authorShim, M.
titleSub-ångström-resolution diffractive imaging
  of single nanocrystals.
journalNature Physics
  volume5, pages129–133
  (year2009).

DeCaro2010
authorDe Caro, L., authorCarlino, E.,
  authorCaputo, G., authorCozzoli, P. D. &
  authorGiannini, C.
titleElectron diffractive imaging of oxygen atoms in
  nanocrystals at sub-angstrom resolution.
journalNat Nano volume5,
  pages360–365 (year2010).
<http://dx.doi.org/10.1038/nnano.2010.55>.

Spence1988
authorSpence, J. C.
titleExperimental high-resolution electron
  microscopy (publisherOxford University Press,
  year1988).

Decaro2013
authorDe Caro, L., authorCarlino, E.,
  authorSiliqi, D. & authorGiannini, C.
titleCoherent diffractive imaging: From nanometric down to
  picometric resolution.
In booktitleHandbook of Coherent-Domain Optical
  Methods, pages291–314 (publisherSpringer,
  year2013).

Fienup1982
authorFienup, J. R.
titlePhase retrieval algorithms: a comparison.
journalAppl. Opt. volume21,
  pages2758–2769 (year1982).
<http://ao.osa.org/abstract.cfm?URI=ao-21-15-2758>.

Abbey2008
authorAbbey, B. et al.
titleKeyhole coherent diffractive imaging.
journalNature Physics
  volume4, pages394–398
  (year2008).

DeCaro2012
authorDe Caro, L., authorCarlino, E.,
  authorVittoria, F. A., authorSiliqi, D. &
  authorGiannini, C.
titleKeyhole electron diffractive imaging (KEDI).
journalActa Crystallographica Section A
  volume68, pages687–702
  (year2012).
<http://dx.doi.org/10.1107/S0108767312031832>.

shechtman2015
authorShechtman, Y. et al.
titlePhase retrieval with application to optical imaging:
  a contemporary overview.
journalIEEE signal processing magazine
  volume32, pages87–109
  (year2015).

Gerchberg1972
authorGerchberg, R. & authorSaxton, W.
titleA practical algorithm for the determination of the
  phase from image and diffraction plane pictures.
journalOptik (Jena) volume35,
  pages237 (year1972).

Marchesini2007
authorMarchesini, S.
titleInvited article: A unified evaluation of iterative
  projection algorithms for phase retrieval.
journalReview of Scientific Instruments
  volume78, pages– (year2007).
<http://scitation.aip.org/content/aip/journal/rsi/78/1/10.1063/1.2403783>.

chen2007
authorChen, C.-C., authorMiao, J.,
  authorWang, C. & authorLee, T.
titleApplication of optimization technique to
  noncrystalline x-ray diffraction microscopy: Guided hybrid input-output
  method.
journalPhysical Review B
  volume76, pages064113
  (year2007).

Goldberg1989
authorGoldberg, D. E.
titleGenetic Algorithms in Search, Optimization and
  Machine Learning (publisherAddison-Wesley Longman Publishing
  Co., Inc., addressBoston, MA, USA, year1989),
  edition1st edn.

thust1997
authorThust, A., authorLentzen, M. &
  authorUrban, K.
titleThe use of stochastic algorithms for phase retrieval
  in high resolution transmission electron microscopy.
journalScanning Microscopy
  volume11, pages437–454
  (year1997).

nicholson1999
authorNicholson, J., authorOmenetto, F.,
  authorFunk, D. & authorTaylor, A.
titleEvolving FROGS: phase retrieval from
  frequency-resolved optical gating measurements by use of genetic algorithms.
journalOptics letters
  volume24, pages490–492
  (year1999).

taylor2006
authorTaylor, J. R., authorKing III, B. A.,
  authorSteincamp, J. & authorRakoczy, J.
titleGenetic algorithm phase retrieval for the systematic
  image-based optical alignment test bed.
journalPublications of the Astronomical Society of
  the Pacific volume118, pages319
  (year2006).

li2011
authorLi, N., authorGao, P., authorLu, Y.,
  authorYu, W. & authorYu, B.
titlePhase retrieval for hard x-ray in-line phase contrast
  imaging based on a parallel hybrid genetic algorithm.
In booktitleComputational Sciences and Optimization
  (CSO), 2011 Fourth International Joint Conference on,
  pages66–70 (organizationIEEE, year2011).

ong2010
authorOng, Y.-S., authorLim, M. H. &
  authorChen, X.
titleResearch frontier-memetic computation–past, present
  & future.
journalIEEE Computational Intelligence Magazine
  volume5, pages24 (year2010).

goldberg1989sizing
authorGoldberg, D. E.
titleSizing populations for serial and parallel genetic
  algorithms.
In booktitleProceedings of the 3rd international
  conference on genetic algorithms, pages70–79
  (organizationMorgan Kaufmann Publishers Inc.,
  year1989).

moscato1989
authorMoscato, P. et al.
titleOn evolution, search, optimization, genetic
  algorithms and martial arts: Towards memetic algorithms.
journalCaltech concurrent computation program, C3P
  Report pages826 (year1989).

renders1996
authorRenders, J.-M. & authorFlasse, S. P.
titleHybrid methods using genetic algorithms for global
  optimization.
journalIEEE Transactions on Systems, Man, and
  Cybernetics, Part B (Cybernetics) volume26,
  pages243–258 (year1996).

el2006
authorEl-Mihoub, T. A., authorHopgood, A. A.,
  authorNolle, L. & authorBattersby, A.
titleHybrid genetic algorithms: A review.
journalEngineering Letters
  volume13, pages124–137
  (year2006).

Neri2012
authorNeri, F. & authorCotta, C.
titleMemetic algorithms and memetic computing
  optimization: A literature review.
journalSwarm and Evolutionary Computation
  volume2, pages1–14 (year2012).

Marchesini2003
authorMarchesini, S. et al.
titleX-ray image reconstruction from a diffraction pattern
  alone.
journalPhys. Rev. B volume68,
  pages140101 (year2003).
<http://link.aps.org/doi/10.1103/PhysRevB.68.140101>.

Smith2007
authorSmith, J. E.
titleCoevolving memetic algorithms: a review and progress
  report.
journalIEEE Transactions on Systems, Man, and
  Cybernetics, Part B (Cybernetics) volume37,
  pages6–17 (year2007).

Millane2015
authorMillane, R. P. & authorArnal, R. D.
titleUniqueness of the macromolecular crystallographic
  phase problem.
journalActa Crystallographica Section A: Foundations
  and Advances volume71 (year2015).
  
SIPI1997
authorWeber, A. G.
titleThe USC-SIPI image database version 5
journalUSC-SIPI Rep. volume315,
pages1–24 (year1997).

Banerjee2015
authorBanerjee, H., authorBanerjee, S.,
  authorRanderia, M. & authorSaha-Dasgupta, T.
titleElectronic structure of oxide interfaces: A
  comparative analysis of GdTiO3/SrTiO3 and LaAlO3/SrTiO3 interfaces.
journalScientific reports
  volume5 (year2015).

Varela2006
authorVarela, M. et al.
titleAtomic scale characterization of complex oxide
  interfaces.
journalJournal of materials science
  volume41, pages4389–4393
  (year2006).

Muller2004
authorMuller, D. A., authorNakagawa, N.,
  authorOhtomo, A., authorGrazul, J. L. &
  authorHwang, H. Y.
titleAtomic-scale imaging of nanoengineered oxygen vacancy
  profiles in SrTiO3.
journalNature volume430,
  pages657–661 (year2004).

Carlino2014
authorCarlino, E.
titleTem for characterization of semiconductor
  nanomaterials.
In booktitleTransmission Electron Microscopy
  Characterization of Nanomaterials, pages89–138
  (publisherSpringer, year2014).

Spence2013
authorSpence, J. C.
titleHigh-resolution electron microscopy
  (publisherOUP Oxford, year2013).

DeCaro2016
authorDe Caro, L., authorScattarella, F. &
  authorCarlino, E.
titleDetermination of the projected atomic potential by
  deconvolution of auto-correlation function of TEM electron nano-diffraction
  patterns.
journalCrystals .
  volume6, pages141
  (year2016).
  
Storn1997
authorStorn, R. & authorPrice, K.
titleDifferential evolution–a simple and efficient
  heuristic for global optimization over continuous spaces.
journalJournal of global optimization
  volume11, pages341–359
  (year1997).

Csiszar1991
authorCsiszar, I.
titleWhy least squares and maximum entropy? an axiomatic
  approach to inference for linear inverse problems.
journalAnn. Statist.
  volume19, pages2032–2066
  (year1991).
<http://dx.doi.org/10.1214/aos/1176348385>.

Carlino2005
authorCarlino, E. & authorGrillo, V.
In booktitleProceedings MCEM VII Portoroze (Si),
  pages159 (year2005).





