Phase Retrieval with Application to Optical Imaging

I Introduction

The problem of phase retrieval, namely – the recovery of a function given the magnitude of its Fourier transform - arises in various fields of science and engineering, including electron microscopy, crystallography, astronomy, and optical imaging. Exploring phase retrieval in optical settings, specifically when the light originates from a laser, is natural, because optical detection devices (e.g., ccd cameras, photosensitive films, the human eye) cannot measure the phase of a light wave. This is because, generally, optical measurement devices that rely on converting photons to electrons (current) do not allow direct recording of the phase: the electromagnetic field oscillates at rates $\sim 10^{15}$ Hz, which no electronic measurement devices can follow. Indeed, optical measurement / detection systems measure the photon flux, which is proportional to the magnitude squared of the field, not the phase. Consequently, measuring the phase of optical waves (electromagnetic fields oscillating at $10^{15}$ Hz and higher) involves additional complexity, typically by interfering it with another (known) field, in the process of holography.

Interestingly, electromagnetic fields do have some other features which make them amenable for algorithmic phase retrieval: their far-field corresponds to the Fourier transform of their near-field. More specifically, given a “mask” that superimposes some structure (an image) on a quasi-monochromatic coherent field at some plane in space, the electromagnetic field structure at a large enough distance from that plane is given by the Fourier transform of the image multiplied by a known quadratic phase factor. Thus, measuring the far-field, magnitude and phase, would facilitate recovery of the optical image (the wavefield). However, as noted above, the optical phase cannot be measured directly by an electronic detector. Here is where algorithmic phase retrieval comes into play, offering a means for recovering the phase given the measurement of the magnitude of the optical far-field and some prior knowledge.

The purpose of this review article is to provide a contemporary review of phase retrieval in optical imaging. It begins with historical background section that also explains the physical setting, followed by a section on the mathematical formulation of the problem. The fourth section discusses existing algorithms, while the fifth section describes various contemporary applications. The last section discusses additional physical settings where algorithmic phase retrieval is important, identifies current challenges and provides a long term vision. This review article provides a contemporary overview of phase retrieval in optical imaging, linking the relevant optical physics to the signal processing methods and algorithms. Our goal is to describe the current state of the art in this area, identify challenges, and suggest vision and areas where signal processing methods can have a large impact on optical imaging and on the world of imaging at large, with applications in a variety of fields ranging from biology and chemistry to physics and engineering.

II Historical Background

Algorithmic phase retrieval offers an alternative means for recovering the phase structure of optical images, without requiring sophisticated measuring setups as in holography. These approaches typically rely on some advanced information in order to facilitate recovery. Back in 1952, Sayre envisioned, in the context of crystallography, that the phase information of a scattered wave may be recovered if the intensity pattern at and between the Bragg peaks of the diffracted wave is finely measured . In crystallography, the material structure under study is periodic (a crystal), hence the far-field information contains naturally strong peaks reflecting the Fourier transform of the periodic information. Measuring the fine features in the Fourier transform enabled the recovery of the phase in some simple cases. Twenty six years later, in 1978, Fienup developed algorithms for retrieving phases of 2D images from their Fourier modulus and constraints such as non-negativity, and a known support of the image (See Fig. 1). In the early eighties, the idea of phase retrieval created a flurry of follow up work, partly because those times signified great hope for realizing an optical computer, of which phase retrieval was supposed to be a key ingredient. However, in the 1980s and 1990s, with the understanding that an optical computer is unrealistic, the interest in algorithmic phase retrieval diminished. Towards the end of the millennium, optical phase retrieval started to come back into contemporary optics research, with the interest arising from a completely different direction: the community of researchers experimenting with X-ray imaging, where new X-ray sources (undulators and synchrotrons) were developed. The wide-spread interest of this field was mainly generated by the first experimental recording and reconstruction of a continuous diffraction pattern (Fourier magnitude squared) of a non-crystalline (non-periodic) test object by Miao and collaborators in 1999 .

The reasons for the revival of optical phase retrieval, in 1999, were actually quite subtle. One goal of optical imaging systems is to increase resolution, that is, to image smaller and smaller features. But, as known since Ernst Abbe’s work in 1873, the highest attainable resolution in diffraction imaging (so-called the diffraction limit) is comparable to the wavelength of the light. For visible light, this diffraction limit corresponds to fraction of microns. Consequently, features on the molecular scale cannot be viewed with visible light in a microscope. One could argue then, why not simply use electromagnetic waves of a much shorter wavelength, say, in the hard X-ray regime, where the wavelength is comparable to atomic resolution? The reason is that lens-like devices and other optical components in this spectral region suffer from very large aberrations and are very difficult to make due to fact that refractive indices of materials in this wavelength regime are close to one. On the other hand, algorithmic phase retrieval is of course not limited by the quality of lenses; however it requires very low noise detectors.

An additional problem is that as resolution is improved (that is, as voxel elements in the recovered image are smaller in size), the number of photons per unit area must obviously increase to provide a reasonable SNR. This means that the required exposure time to obtain a given signal level must increase as $(1/d)^{4}$, with $d$ being the resolution length, assumed to be larger than atomic scales . This, in turn, creates another problem: X-ray photons are highly energetic. The atomic cross section for photoabsorption is usually much higher than for elastic scattering, meaning that for every photon that contributes to the diffraction pattern (the measured Fourier magnitude), a considerable greater number of photons are absorbed by the sample. This energy dissipates in the sample first by photoionisation and the breakage of bonds, followed by a cascade of collisional ionisation by free electrons and, at longer timescales, a destruction of the sample due to radiolysis, heating, or even ablation of the sample. Such radiation damage hinders the ability to recover the structure of molecules: the measured far-field intensity (Fourier magnitude) also reflects the structural damages, rather than providing information about the true molecular structure . A solution to this problem was suggested by Solem and Chapline in the 1980’s. They proposed to record images (or holograms in their case) with pulses that are shorter than the timescale for the X-ray damage to manifest itself at a particular resolution. They predicted that picosecond pulses would be required to image at nanometer length scales . Towards the late nineties, with the growing promise in constructing X-ray lasers that generate ultrashort pulses on the femtosecond scale, it was suggested that such pulses could even outrun damage processes at atomic length scales. However, forming a direct image in this way would still require high quality optical components (lenses, mirrors) in the X-ray regime, which do not currently exist. This is because creating lenses for the hard X-ray wavelength regime requires fabrication at picometer resolution, much smaller even than the Bohr radius of atoms. Likewise, while mirrors for X-rays do exist, their best resolution is on the scale of many nanometers, much larger than the features one would want to resolve in imaging of molecules, for example.

The difficulties outlined above in direct X-ray imaging leave no choice but to use alternative methods to recover the structure of nanometric samples. Here is where phase retrieval can make its highest impact. Placing an area detector far enough from the sample to record the far-field diffraction intensity (which is approximately proportional to the squared magnitude of the Fourier transform of the image, if the coherence length of the X-ray wave is larger than the sample size ), together with appropriate constraints on the support of the sample, enable the recovery of the image at nanometric resolution. Indeed, the phase information has been shown numerically and experimentally to be retrieved in this fashion in various examples . The combination of X-ray diffraction, oversampling and phase retrieval has launched the currently very active field called Coherent Diffraction Imaging or CDI . In CDI, an object is illuminated with a coherent wave, and the far field diffraction intensity pattern (corresponding to the Fourier magnitude of the object) is measured. The problem then is to recover the object from the measured far-field intensity (See box on Coherent Diffractive Imaging and Fig. 2 within). Since its first experimental demonstration, CDI has been applied to image a wide range of samples using synchrotron radiation , X-ray free electron lasers (XFELs) , high harmonic generation , soft X-ray laser , optical laser , and electrons . Several readable reviews on the development and implementation of phase-retrieval algorithms for the specific application of CDI were written by Marchesini , Thibault and Elser and Nugent . Presently, one of the most challenging problems in CDI is towards 3D structural determination of large protein molecules . There has been ongoing progress towards this goal during the past decade. In 2006, Chapman et al., demonstrated the CDI of a test sample using intense ultra-short single pulse from free electron laser, relying on recording a diffraction pattern before the sample was destroyed . Recently, the technique was implemented for high-resolution imaging of isolated sub-micron objects, such as herpesvirus , mimivirus and aerosol particles such as soot .

From a theoretical and algorithm perspective, phase retrieval is a difficult problem, in many cases lacking a unique solution. Furthermore, even with the existence of a unique solution, there is not necessarily a guarantee that it can be found algorithmically. Nevertheless, as reasoned above, phase retrieval algorithms and applications have benefited from a surge of research in recent years, in large part due to various new imaging techniques in optics. This trend has begun impacting the signal processing community as well – the past few years have witnessed growing interest within this community in developing new approaches to phase retrieval by using tools of modern optimization theory . More recent work has begun exploring connections between phase retrieval and structure-based information processing . For example, it has been shown that, by exploiting the sparsity of many optical images, one can develop powerful phase retrieval methods that allow for increased resolution considerably beyond Abbe’s diffraction limit, resolving features smaller than $1/5$ of the wavelength . The relationship between the fields of sparsity and optical imaging has led to an important generalization of the basic principles of sparsity-based reconstruction to nonlinear measurement systems . Here too, optics played an important role in signal processing: since the phase retrieval problem is inherently mathematically nonlinear (i.e., the sought signal is related to the measurements nonlinearly), employing sparsity-based concepts in phase retrieval required genuine modifications to the linear sparsity-based algorithms known from the field of compressed sensing . We believe that this field will grow steadily in the next few years, with rapid development of coherent X-ray sources worldwide and more researchers contributing to the theory, algorithms and practice of nonlinear sparse recovery.

III Mathematical Formulation

Recovery of $\mathbf{x}$ from measurement of ${\bf X}$ can be achieved by simply applying the inverse-DFT operator. Writing $X[k]=|X[k]|\cdot e^{j\phi[k]}$, the Fourier phase retrieval problem is to recover ${\bf x}$ when only the magnitude of $X$ is measured, i.e. to recover $x[n]$ given $|X[k]|$. Since the DFT operator is bijective, this is equivalent to recovering the phase of $X[k]$, namely, $\phi[k]$ - hence the name phase retrieval. Denote by $\mathbf{\hat{x}}$ the vector $\mathbf{x}$ after padding with $N-1$ zeros. The autocorrelation sequence of $\mathbf{\hat{x}}$ is then defined as:

It is well known that the DFT of $g[m]$, denoted by $G[k]$, satisfies $G[k]=|X[k]|^{2}$. Thus, the problem of recovering a signal from its Fourier magnitude is equivalent to the problem of recovering a signal from the autocorrelation sequence of its oversampled version.

Continuous phase retrieval can be defined similarly to its discrete counterpart, as the recovery of a 1D signal $f(x)$ from its continuous Fourier magnitude:

The actual objects of interest, electromagnetic fields, are usually described by continuous functions. However, since the data acquisition is digitized (by CCD camera and alike), and the processing is done digitally, we shall mostly treat here the discrete case.

The Fourier phase retrieval problem is as a special case of the more general phase retrieval problem, where we are given measurements:

with $\mathbf{a}_{k}$ denoting the measurement vectors. In discrete 1D Fourier phase retrieval the measurement vectors are given by $\mathbf{a}_{k}[n]=e^{-j2\pi\frac{kn}{M}}$. For mathematical analysis, it is often easier to treat the case where the measurements are random (i.e. $\mathbf{a}_{k}$ are random vectors), as this allows uniqueness guarantees that are otherwise hard to obtain . Nevertheless, more structured measurements have also been investigated .

Before proceeding to the mathematical methodology, it is important to highlight the significance of knowing the Fourier phase. In fact, it is well known that knowledge of the Fourier phase is crucial in recovering an object from its Fourier transform . Many times the Fourier phase contains more information than the Fourier magnitude, as can be seen in the synthetic example shown in Fig. 3. The figure shows the result of the following numerical experiment: Two images (Cameraman and Lenna) are Fourier transformed. The phases of their transforms are swapped, and subsequently they are inverse Fourier transformed. It is evident, for this quite arbitrary example, that the Fourier phase contains a significant amount of information about the images. In crystallography, this phenomenon is the source of genuine concern of phase bias of molecular models (such as used in molecular replacement) in refined structures.

In the remainder of this section we discuss uniqueness of the phase retrieval problem, i.e. under what conditions is the solution to the phase problem unique? It is worth noting that, while the discussion of theoretical uniqueness guarantees is important and interesting, the lack of such guarantees does not prevent practical applications from producing excellent reconstruction results in many settings.

III-B Uniqueness

The recovery of a signal from its Fourier magnitude alone, in general, does not yield a unique solution. This section will review the main existing theoretical results regarding phase-retrieval uniqueness.

First, there are so called trivial ambiguities that are always present. The following three transformations (or any combination of them) conserve Fourier magnitude:

1. Global phase shift: $x[n]\Rightarrow x[n]\cdot e^{j\phi_{0}}$

2. Conjugate inversion: $x[n]\Rightarrow\overline{x[-n]}$

3. Spatial shift: $x[n]\Rightarrow x[n+n_{0}]$.

Second, there are non-trivial ambiguities, the situation of which varies for different problem-dimensions. In the 1D problem there is no uniqueness – i.e. there are multiple 1D signals with the same Fourier magnitude. Even if the support of the signal is bounded within a known range, uniqueness does not exist . Any pair of 1D signals having the same autocorrelation function yields the same Fourier magnitude, as the two are connected by a Fourier transform. Consider for example the two vectors $\mathbf{u}=[1\,\,0\,\,-2\,\,0\,\,-2]^{T}$ and $\mathbf{v=}[(1-\sqrt{3})\,\,0\,\,1\,\,0\,\,(1+\sqrt{3})]^{T}$. Both of these vectors have the same support, and yield the same autocorrelation function $g[m]=$. Therefore, they are indistinguishable by their Fourier magnitude, even though they are not trivially equivalent.

For higher dimensions (2D and above), Bruck and Sodin , Hayes , and Bates have shown that, with the exception of a set of signals of measure zero, a real $d\geq 2$ dimensional signal with support $\mathbf{N}=[N_{1}\dots N_{d}]$, namely $x(n_{1},\dots,n_{d})=0$ whenever $n_{k}<0$ or $n_{k}\geq N_{k}$ for $k=1,\dots,d$ is uniquely specified by the magnitude of its continuous Fourier transform, up to the trivial ambiguities mentioned above. Furthermore, the magnitude of the oversampled $\mathbf{M}$ point DFT sequence of the signal, with $\mathbf{M\geq}2\mathbf{N}-1$ (where the inequality holds in every dimension), is sufficient to guarantee uniqueness. The problematic set of signals that are not uniquely defined by their Fourier magnitudes are those having a reducible Z transform: denoting the $d$ dimensional $Z$ transform of $x$ by $X(z_{1},\dots,z_{d})=\sum_{n_{1}}\cdots\sum_{n_{d}}x(n_{1},\dots,n_{d})z_{1}^{-n_{1}}\cdots z_{d}^{-n_{d}}$ , $X(\mathbf{z})$ is said to be reducible if it can be written as $X(\mathbf{z})=X_{1}(\mathbf{z})X_{2}(\mathbf{z})$, where $X_{1}(\mathbf{z})$ and $X_{2}(\mathbf{z})$ are both polynomials in $\mathbf{z}$ with degree $p>0$. It is important to note that in practice, for typical images, a number of samples smaller than $2\mathbf{N}-1$ is many times sufficient (even $\mathbf{M=N}$ can work ), however the exact guarantees relating the number of samples to the type of images remains an open question.

III-B2 General measurements

The results reviewed in this section are summarized in Table I. In addition, there is a large amount of work on phase retrieval uniqueness under different conditions, e.g. when the phase is known only approximately , or from redundant masked Fourier measurements .

IV Algorithms

Despite the uniqueness guarantees, no known general solution method exists to actually find the unknown signal from its Fourier magnitude given the other constraints. Over the years, several approaches have been suggested for solving the phase retrieval problem, with the popular ones being alternating projection algorithms . In addition, methods were suggested attempting to solve phase retrieval problems by using exposures with different masks , or images obtained at different propagation planes . Another method to obtain additional information is scanning CDI, (also termed ptychography) , which uses several different illumination patterns to obtain coherent diffraction images.

In this section we survey existing phase retrieval algorithms, including general algorithms (sub-Section IV-A), and sparsity based algorithms, i.e. algorithms exploiting prior knowledge in the form of signal-sparsity (sub-Section IV-B).

The general phase retrieval problem we wish to solve can be formulated as the following least squares problem, or empirical risk minimization:

with $\mathbf{y}$ being the measurements and $\mathbf{a}_{k}$ being the measurement vectors defined in (4). In general we can replace the square in the objective by any power $p$. Unfortunately, this is a non-convex problem, and it is not clear how to find a global minimum even if one exists. In this section we describe several approaches that have been suggested to deal with this problem, and types of prior information that can be incorporated into these methods in order to increase the probability of convergence to the true solution.

The most popular class of phase retrieval methods is based on alternate projections. These methods were pioneered by the work of Gerchberg and Saxton , dealing with the closely related problem of recovering a complex image from magnitude measurements at two different planes - the real (imaging) plane and Fourier (diffraction) plane. The original Gerchberg-Saxton (GS) algorithm consists of iteratively imposing the real-plane and Fourier-plane constraints, namely, the measured real-space magnitude $|x[n]|$ and Fourier magnitude $|X[k]|$, as illustrated in Fig. 4a. The GS algorithm is described in Algorithm 1. The recovery error, defined as E_{i}=\sum_{k}\Big{|}|Z_{i}[k]|-|X[k]|\Big{|}^{2} is easily shown to be monotonically decreasing with $i$ . Despite this fact, recovery to the true solution is not guaranteed, as the algorithm can converge to a local minimum.

Extending the GS projection ideas further, Fienup in 1978 suggested a modified version of the GS algorithm, in which the real-space magnitude constraints may be replaced by other types of constraints, in addition to consistency with the measured Fourier magnitude. The real-space constraints may be for example non-negativity, a known signal support, namely $x[i]=0$ for all $i>N_{0}$, where $N_{0}$ is known (or approximately known), or both. The basic framework of the Fienup methods is similar to GS - in fact, the first three steps are identical. Step 4), however, replaces imposing the real-space magnitude constraint by applying a correction to the real-space estimate. Some possible variants to this step were also suggested . Here we describe the one most commonly used, referred to as the hybrid-input output (HIO) method, which consists of the following correction step:

Obtain $z_{i+1}[n]$ by applying a correction to the real-space image estimate:

with $\beta$ being a small parameter, and $\gamma$ being the set of indices for which $z_{i}^{\prime}[n]$ violates the real-space constraints.

The real-space constraint violation may be a support violation (signal is nonzero where it should be zero), or a non-negativity violation.

The Fienup algorithm is represented schematically in Fig. 4b. There is no proof that the HIO algorithm converges. It is also known to be sensitive to the accuracy of the prior information (e.g. the real-space support needs to be tightly known, especially in the complex signal case ). Nonetheless, in practice, the simple HIO based methods are commonly used in optical phase retrieval applications such as CDI . Other variants of the correction step include the Input-Output method, and the Output-Output method , corresponding respectively to

An important feature of the HIO algorithm is its empirical ability to avoid local minima and converge to a global minimum for noise-free oversampled diffraction patterns. However, when there is high noise present in the diffraction intensity, HIO suffers from several limitations. First, the algorithm sometimes becomes stagnant and fails to converge to a global minimum. Second, a support has to be pre-defined. Third, the image oscillates as a function of the iteration. Over the years, various algorithms have been developed to overcome these limitations, including the combination of HIO and the error-reduction (ER) algorithm , difference map , hybrid projection reflection , guided hybrid input-output (GHIO) , relaxed averaged alternating reflectors (RAAR) , noise robust (NR)-HIO , and oversampling smoothness (OSS) . An analysis of iterative phase retrieval algorithms from a convex optimization perspective can be found in .

As an exapmle, the recently proposed OSS algorithm exhibits improved performance over HIO and its variants in many settings. OSS is based on Fienup iterations, with an added smoothing Gaussian filter applied to the off-support region in the real space object in each iteration. Namely, the fourth step in HIO is replaced by:

where $W[k]$ is a Gaussian function, with its variance decreasing with iterations. A quantitative comparison for a specific example between OSS and HIO can be found in Section V-A. For a comparison and numerical investigation of several alternate projection algorithms see for example .

As performance of the Fienup methods is dependent on the initial points, it is possible and recommended to try several initializations. In , the authors consider a clever method for initial point selection, and show that for the random Gaussian measurement case, the resulting iterations yield a solution arbitrarily close to the true vector.

IV-A2 Semi-Definite Programming (SDP) based algorithms

An alternative recently developed to solve the phase retrieval problem is based on semidefinite relaxation . The method relies on the observation that (4) describes a set of quadratic equations, which can be re-written as linear equations in a higher dimension. Specifically, define the $N\times N$ matrix ${\bf X}={\bf x}{\bf x}^{*}$. The measurements (4) are then linear in ${\bf X}$:

where ${\bf A}_{k}=\mathbf{a}_{k}\mathbf{a}_{k}^{*}$. Our problem is then to find a matrix ${\bf X}={\bf x}{\bf x}^{*}$ that satisfies (9). The constraint ${\bf X}={\bf x}{\bf x}^{*}$ is equivalent to the requirement that ${\bf X}$ has rank one, and is positive semi-definite, which we denote by ${\bf X}\succeq 0$ . Therefore, finding a vector ${\bf x}$ satisfying (4) can be formulated as:

Problem (10) is equivalent to the following rank minimization problem:

Unfortunately, rank minimization is a hard combinatorial problem. However, since the constraints in (11) are convex (in fact linear), one might try to relax the minimum rank objective, for example by replacing it with minimization of $\textnormal{Tr}({\bf X})$. This approach is referred to as PhaseLift . Alternatively, one may use the log-det reweighted rank minimization heuristic suggested in , which is the approach followed in . In it is shown that PhaseLift yields the true vector ${\bf x}$ with large probability, when the measurements are random Gaussian and $M\sim O(N\log N)$.

The SDP approach requires matrix lifting, namely, replacing the sought vector with a higher dimensional matrix, followed by solving a higher dimensional problem. It is therefore, in principle, more computationally demanding than the alternating projection approaches, or greedy methods, which will be discussed in the next section. In addition, in general there is neither a guarantee that the rank minimization process will yield a rank-one matrix, nor that the true solution is found, even if there is a unique solution.

IV-B Sparsity based algorithms

A specific kind of prior knowledge that can be incorporated into the phase retrieval problem to help regularize it, is the fact that the sought real-space object is sparse in some known representation (See sparse linear problems box). This means that the object $\mathbf{x}$ can be written as:

with $\mathbf{\boldsymbol{\Psi}}$ being a representation matrix (the sparsity basis), and $\boldsymbol{\alpha}$ being a sparse vector, i.e. a vector containing a small number of nonzero coefficients. The simplest example is when the object is composed of a small number of point sources (in which case $\boldsymbol{\Psi}$ is the identity matrix). Armed with such prior knowledge, one can hope to improve the performance of phase retrieval algorithms, by limiting the search for the true vector to the set of sparse vectors. There are several different ways that such knowledge can be incorporated, which are described in this subsection.

The Fienup algorithm, described in Section IV-A1, allows in principle incorporation of various types of general knowledge about the object, including sparsity . The method in , for example, is based on the Fienup iterations, with the first three steps remaining unchanged. Step 4, is replaced by projection and thresholding. Assuming an invertible $\boldsymbol{\Psi}$ and a $k$ sparse vector $\boldsymbol{\alpha}$ such that ${\bf x}=\boldsymbol{\Psi\alpha}$:

Obtain $z_{i+1}[n]$ by projecting $z_{i}^{\prime}[n]$ onto $\boldsymbol{\Psi}$, thresholding, and projecting back. Namely:

Calculate $\boldsymbol{\alpha}{}_{i}=\boldsymbol{\Psi}^{-1}\mathbf{z}_{i}^{\prime}$

Keep only the $k$ largest elements of $\mathbf{|a}_{i}|$, setting the rest to zero.

Set $\mathbf{z}_{i+1}=\boldsymbol{\Psi}\boldsymbol{\alpha}_{i}$.

Similarly to the GS method, the error here can be shown to be nonincreasing, so that convergence to a local minimum is guaranteed .

Note, that while this method is suggested in for an orthonormal basis $\boldsymbol{\Psi}$, it can be easily modified to accommodate a non-invertible $\boldsymbol{\Psi}$. This can be done by replacing parts (a)+(b) with finding a sparse solution $\mathbf{a}_{i}$ to $\mathbf{z}_{i}=\boldsymbol{\Psi}\mathbf{a}_{i}$, using any sparse solution method .

SDP based methods with sparsity prior

SDP based methods can also be modified to account for prior knowledge of signal sparsity. The incorporation of sparsity can be performed in several different ways. The first work to suggest sparsity-based SDP phase-retrieval came from the domain of optics, and dealt with partially spatially-incoherent illumination . This work actually considered a theoretical problem of greater complexity, combining phase retrieval with sub-wavelength imaging. Experimental results on sub-wavelength CDI can be found in , where the sought signal is an actual optical image with subwavelength features, and the measured data corresponds to the Fourier magnitude sampled by a camera at the focal plane of a microscope lens.

The method suggested in , dubbed QCS for Quadratic Compressed Sensing, is based on adding sparsity constraints to the rank minimization problem (11). When ${\bf x}$ is sparse, the result of the outer product ${\bf X}={\bf x}{\bf x}^{*}$ is a sparse matrix as well, as shown in Fig. 5. Therefore, one strategy might be to minimize the $l_{1}$ norm of the matrix ${\bf X}.$ Alternatively, it is possible to exploit further the structure of ${\bf X},$ by noticing that the number of rows in ${\bf X}$ with a nonzero norm is equal to the number of non-zero values in ${\bf x}$. This means that sparsity of $\mathbf{x}$ also implies a small number of non-zero rows in ${\bf X}$. Consider the vector $\mathbf{p}$ containing the $l_{2}$ norm of the rows of ${\bf X}$, i.e. p_{j}=\big{(}\sum_{k}|X_{jk}|^{2}\big{)}^{\frac{1}{2}} (note that the $l_{2}$ norm can be replaced by any other norm). Since $\mathbf{p}$ should be sparse, one might try to impose a low $l_{1}$ norm on $\mathbf{p}$, in the spirit of $l_{1}$ minimization for the sparse linear problem. This yields the constraint ||\mathbf{p}||_{1}=\sum_{j}|p_{j}|=\sum_{j}\big{(}\sum_{k}|X_{jk}|^{2}\big{)}^{\frac{1}{2}}\leq\eta, corresponding exactly to a low mixed $l_{1-2}$ norm constraint on ${\bf X}$ . The problem to solve, as cast in is therefore :

where $\epsilon$ is a noise parameter, and $\eta$ is a sparsity parameter, enforcing row sparsity of $\mathbf{X}$.

Since finding a rank 1 matrix $\mathbf{X}$ satisfying the constraints is NP hard, the solution to (19) is approximated in () using the iterative log-det heuristic proposed in , with an additional thresholding step added at each iteration, to further induce signal sparsity. Once a low rank matrix $\hat{\mathbf{X}}$ that is consistent with the measurements and the sparse priors is found, the sought vector $\mathbf{x}$ is estimated by taking the best rank 1 approximation of $\hat{{\bf X}}$ using the singular value decomposition (SVD): Decomposing $\hat{{\bf X}}$ into $\hat{{\bf X}}=\mathbf{USV^{T}}$, the rank-1 approximation of $\hat{{\bf X}}$ is taken as $\mathbf{\hat{X}_{1}}=S_{11}\mathbf{U_{1}}\mathbf{U_{1}}^{*}$, where $S_{11}$ represents the largest singular value, and $\mathbf{U_{1}}$ is the corresponding column of $\mathbf{U}$.

Similar ideas that add sparse priors to SDP methods have been later suggested in . In , the rank minimization objective is relaxed to a convex trace minimization, with an additional $l_{1}$ regularization term to induce sparsity. This formulation yields:

The solution of (20) is shown to be unique in the noiseless case ($\epsilon=0$), under the following condition: ||\bar{{\bf X}}||_{0}\leq\frac{1}{2}\big{(}1+\frac{1}{\mu}), where $\bar{{\bf X}}=\bar{{\bf x}}\bar{{\bf x}}^{*}$, with $\bar{{\bf x}}$ being the true solution to (4). The mutual coherence $\mu$ is defined by $\mu=\max_{i,j}\frac{<\mathbf{B}_{i},\mathbf{B}_{j}>}{||\mathbf{B}_{i}||||\mathbf{B}_{j}||}$ , with $\mathbf{B}$ being the matrix satisfying ${\bf y}=\mathbf{BX^{S}}$, where $\mathbf{X^{S}}$ is the vector obtained from stacking the columns of $\mathbf{X}$. The same work also relates other recovery guarantees to the RIP criterion.

In it is shown that for $\mathbf{a}_{i}$ that are independent, zero-mean normal vectors, on the order of $k^{2}\log n$ measurements are sufficient to recover a $k$-sparse input from measurements of the form (4), using SDP relaxation. In , an algorithm is suggested to solve the sparse 1D Fourier phase retrieval problem based on a two-step process, each step cast separately as an SDP problem: First, the support of ${\bf x}$ is determined from its autocorrelation sequence, and then ${\bf x}$ is determined, given the support. This algorithm is shown experimentally to recover $k$ sparse signals from $O(k^{2})$ measurements.

Greedy methods with sparsity prior

Since the matrix lifting algorithms involve a dimension increase, they are not ideally suited for large vectors, where computational cost can become significant. In addition, they are in general not guaranteed to converge to a solution. An alternative to the SDP algorithms is posed by sparsity based greedy algorithms . One algorithm, that is both fast and accurate, is a greedy method named GESPAR (for GrEedy Sparse PhAse Retrieval) . GESPAR attempts to solve the least squares sparse quadratic problem (5). Namely, it seeks a $k$ sparse vector $\mathbf{x}$ consistent with the quadratic measurements $\mathbf{y}$. GESPAR is a fast local search method, based on iteratively updating signal support, and seeking a vector that corresponds to the measurements, under the current support constraint. A local-search method is repeatedly invoked, beginning with an initial random support set. Then, at each iteration a swap is performed between a support and an off-support index. Only two elements are changed in the swap (one in the support and one in the off-support), following the so-called 2-opt method . Given the support of the signal, the phase-retrieval problem is then treated as a non-convex optimization problem, approximated using the damped Gauss Newton method . See Algorithm 2 for a detailed description of the algorithm.

GESPAR has been shown to yield fast and accurate recovery results (see Sparse phase retrieval algorithms box and Fig. 7 within), and has been used in several phase-retrieval optics applications - including CDI of 1D objects , efficient CDI of sparsely temporally varying objects , and phase retrieval via waveguide arrays . A similar method has been used to solve the combined phase-retrieval and sub-wavelength imaging problem (See sub-Section V-D).

V Applications in Lensless Imaging

In this section, we present several CDI applications with connection to the phase retrieval algorithms described in previous sections. The concept of phase retrieval in optical imaging arises from the attempt to recover images from experimental measurements. To this end, it is essential to emphasize that compared to numerical simulations or signal processing of digital data, phase retrieval of experimentally obtained patterns has several additional challenges. First, the far-field intensity distribution (Fourier magnitude) is corrupted by various types of noise, such as Poisson noise, detector read-out noise, and unwanted parasitic scattering from the optics components in the system. Second, in single-shot experiments, the measured far-field intensity distribution is usually incomplete, including a missing center (i.e. the very low spatial frequency information cannot be directly recorded by a detector) . Third, when the far-field intensity distribution is measured by a detector, each pixel integrates the total number of photons within the solid angle subtended by the pixel, which is not exactly equivalent to uniform sampling of the diffraction signal . Additionally, many experiments are carried out using incoherent (but bright) sources. Coherence is achieved by propagating a long distance from the source, but often the experiment is constrained to be carried out with a partially-incoherent beam . All of these issues add complications to algorithmic phase retrieval. However, notwithstanding these challenges, successful phase retrieval of experimental data in optical imaging has been widely achieved . Below, we show several examples.

Quantitative comparisons between the OSS, HIO, ER-HIO and NR-HIO algorithms have been performed using both simulated and experimental data . Figure 8 shows a noise-free oversampled diffraction pattern (Fourier magnitude squared) calculated from a simulated biological vesicle (Fig. 8c). High Poisson noise was then added to the diffraction intensity (Fig. 8b). Figures 8d-g show the final reconstructions by HIO, ER-HIO, NR-HIO, and OSS, respectively. Visually, OSS produced the most faithful reconstruction among the four algorithms (insets in Fig. 8d-g). The recovery error was quantified using consistency with the measurements:

where $z_{r}[n]$ is the final reconstruction and $z_{m}[n]$ is the model structure. The value for $E$ of the HIO, ER-HIO, NR-HIO, and OSS reconstructions is $0.28,0.24,0.16$ and $0.07$, respectively.

Next, the four algorithms were compared using an experimental diffraction pattern measured from a Schizosaccharomyces pombe yeast spore cell . The experiment was conducted on an undulator beamline at a 3rd generation synchrotron radiation (Spring-8) in Japan. A coherent wave of 5 keV X-rays was incident on a fixed, unstained S. pombe yeast spore. An oversampled X-ray diffraction pattern was acquired by a charge-coupled device detector. Figure 9a shows the experimental diffraction pattern in which the centro-square represents the missing low spatial resolution data . By using a loose support, phase retrieval was performed on the measured data with the HIO, ER-HIO, NR-HIO, and OSS algorithms. For each algorithm, five independent trials were conducted, each consisting of 100 independent runs with different random initial phase sets. In each trial, the reconstruction with the smallest error metric $R_{F}$ was chosen as a final image, where $R_{F}$ is defined as:

Here $|Z_{e}[k]|$ is the measured Fourier magnitude, $|Z_{m}[k]|$ is the recovered Fourier magnitude, and $\zeta$ is a scaling factor.

For each algorithm, the mean and average of the five final images were used to quantify the reconstruction. Figures 9c-j show the average and variance of five final images obtained by HIO, ER-HIO , NR-HIO , and OSS , respectively. The average $R_{F}$ and the consistency of five independent trials are shown in Fig. 9b. Both visual inspection and quantitative results indicate that OSS produced the most consistent reconstructions among all four algorithms.

V-B X-ray free electron laser CDI

The majority of imaging experiments at X-ray free-electron laser sources utilize the method of CDI. The lensless nature is actually an advantage when dealing with extremely intense and destructive pulses, where one can only carry out a single pulse measurement with each object (say, a molecule) before the object disintegrates. In such cases, often one cannot use any optical components at all, because any component (e.g., a lens) would be severely damaged by the extremely high flux of (X-ray) photons. CDI solves these problems: it works without the need for optical components. In this vein, CDI also facilitates reliable imaging of moving objects. Indeed, in many experiments the objects move (flow) across the X-ray beam, for example, when the X-ray laser beam hits a focused aerosol beam or nano-particles in a liquid jet. In such an experiment, the particle density is usually adjusted so that the X-ray laser pulse is more likely to hit a single particle than several. A particle is hit by chance by a pulse, but this is not known until the diffraction pattern is read out from the detector, which is done on every pulse. The stream of data is then analyzed and sorted to give the single-particle hits, which contain the meaningful measured data, while all other data is ignored.

There are two generic classes of these “single particle” CDI experiments: imaging of reproducible particles, and imaging of unique particles. The first category includes particles such as viruses. Assuming that these particles are not aligned in the same direction, the collected data represents diffraction patterns of a common object, but in random orientations. If the orientations can be determined then the full three-dimensional Fourier magnitude of the object can be determined, which in turn could be phased to give a 3D image. A proof of concept of this experiment was carried out by Loh et al .

An example of the second class of flash diffractive imaging is imaging airborne soot particles in flight in an aerosol beam . Several diffraction patterns of soot particles and clusters of polystyrene spheres (as test objects) are shown in Fig. 10, along with the 2D reconstructions of the objects. The experiments were carried out at the LCLS, using the CFEL-ASG Multi-Purpose (CAMP) instrument at the Atomic, Molecular and Optical Science beam line . Pulses of about $10^{12}$ photons of 1.0 nm wavelength were focused to an area of 10 $\mu$$m^{2}$. The X-ray detectors (pnCCD panels) were placed to give a maximum full-period resolution of 13 nm at their center edges.

In these experiments, the phase retrieval of the patterns was carried out using the Relaxed Averaged Alternating Reflections (RAAR) algorithm, and using the Shrinkwrap procedure , which determines and iteratively updates the support constraint used. The objects were such that it was possible to apply an additional constraint that the image is real valued. Strikingly, the X-ray coherent diffraction patterns have very high contrast. The intensity minima are close to zero. This has an enormous effect on the ability to recover the phase of these pattern reliably. This reliability is quantified in the phase-retrieval transfer function (PRTF) , which compares the magnitude of the complex-valued average of patterns phased with different starting guesses to the square root of the measured diffraction pattern. If, at a particular pixel of the diffraction pattern, the phases are consistently reconstructed, then the sum over $N$ patterns will give a magnitude $N$ times higher than the measured magnitude, and so the PRTF will be unity. If the phases are random, then this sum will tend to zero. For patterns generated with X-Ray free electron lasers, this function is often close to unity and is lower primarily in areas where the signal to noise is low. This limited signal is what ultimately limits the resolution; an estimate of the achieved resolution is given by the white dotted circle on each pattern in Fig. 10. The reconstructed images are sums of ten independent reconstructions. These complex-valued sums have the nice property that the their Fourier spectrum is effectively modulated by the PRTF and hence any artifacts due to noise (or even due to forcibly truncating the data to a lower resolution) is unlikely to show up in the recovered image.

V-C Tabletop short wavelength CDI

To-date, most CDI experiments are carried out in 3rd generation synchrotron and X-ray free electron lasers. However, limited access and experimental time hinder the development and applications of CDI using these methods. Thus, over the past several years, CDI microscopes that are based on tabletop sources of coherent extreme UV and soft X-rays are also being developed . Figure 11 shows the first tabletop CDI experiment with extreme UV wavelength.

V-D Sub-wavelength CDI using sparsity

Prior knowledge of object-sparsity can help regularize the phase-retrieval problem, as well as compensate for loss of other kinds of information - in this example - the loss of high spatial frequencies. As described before, when an object is illuminated by coherent light of wavelength $\lambda$, the far-field intensity pattern is proportional to the magnitude of the object’s Fourier transform. In addition, features in the object that are smaller than $\sim\lambda/2$ are smeared due to the diffraction limit. Consequently, the intensity measured in the far field corresponds to $\mathbf{y}\propto|\mathbf{LF}\mathbf{x}|^{2}$ where $\mathbf{L}$ represents a low pass filter at cutoff frequency $\nu_{c}=1/\lambda$, $\mathbf{F}$ represents the Fourier transform, and $|\cdot|^{2}$ stands for elementwise squared absolute value.

Figure 12 (adapted from ), shows the recovery of a sparse object containing sub-wavelength features ($100nm$ holes illuminated by a $\lambda=532nm$ laser) from its experimentally measured low pass filtered Fourier magnitude. The prior knowledge used for recovery is that the object is comprised of a small number of $100nm$ diameter circles on a grid, illuminated by a plane wave. The exact number, locations, and amplitudes of the circles are not known a-priori. The recovery is performed using a greedy algorithm that iteratively updates the support of the object, finds a local minimum and removes the weakest circle, until convergence .

Another type of information loss in CDI, for which the prior knowledge of object sparsity can be helpful, is low signal to noise ratio. In non-destructive X-ray CDI measurements, it is not uncommon for signal acquisition time to be on the order of tens of seconds , in order to acheive sufficiently high SNR. This poses a severe limitation on the temporal resolution attainable with such measurements, restricting the types of dynamical phenomena accessible by X-ray CDI. Expoilting sparsity in the change that an object undergoes between subsequent CDI frames has been recently suggested as a means to overcome high noise values, and consequently significantly decrease acquisition time . In other words, the fact that an object is sparsely varying, can be used as prior information to effectively denoise sequential Fourier magnitude measurements. In , CDI of a sparsely varying object is formulated as a sparse quadratic optimization problem, and solved using GESPAR . Numerical simulations suggest a dramatic potential improvement in temporal resolution: In an example consisting of a $51\times 51$ pixel object with 5 randomly varying pixels between frames, an improvement of 2 orders of magnitude in acquisition time is possible .

VI Other Physical Settings, bottlenecks, and Vision

This review article focused mainly on the simplest physical setting for phase retrieval in optical imaging (Fig. 2), CDI: where an unknown 2D optical image is recovered algorithmically from a single measurement of its far-field intensity pattern, given a known image support (or other prior information). In terms of signal processing, this problem corresponds to recovering a 2D object from measurements of its Fourier magnitude. However, the issue of phase retrieval in optical imaging, and in a more general sense – in optics, is far broader, and includes other physical settings which naturally translate into signal processing problems different than the standard phase retrieval problem. This section provides a short overview of those physical settings, defines the various problems in terms of signal processing, and provides some key references. We conclude with a discussion on the main challenges and bottlenecks of phase retrieval in optical imaging, followed by an outlook for the upcoming years and long term vision.

The simplest optical phase retrieval problem assumes that the measured data corresponds to the Fourier magnitude. In optical settings, this means that the measurements are taken in the Fourier domain of the sought image, which physically means performing the measurements at a plane sufficiently far away from the image plane (the so-called far-field or the Fraunhofer regime), or at the focal plane of an ideal lens . In reality, however, the measurements can be taken at any plane between the image plane and the far field, which would yield intensity patterns that are very different than the Fourier magnitude of the sought signal. This of course implies that new (or revised) algorithms - beyond those described in previous sections - must be used, which naturally raises issues of conditions for uniqueness and convergence. At the same time, these measurements have some interesting advantages, which can be used wisely to improve the performance of phase retrieval. Let us begin by describing the relevant physical settings.

As stated earlier in this paper, the optical Fourier plane corresponds to a plane sufficiently far away from where the object (the sought signal) is positioned. Far away here means asymptotically at infinite distance from the object plane, or at the focal plane of a lens. However, the entire propagation-evolution of electromagnetic waves from any plane to any other plane (not only from the near field to the far field) is known: it is fully described by Maxwell’s equations. As such, one can formulate the problem through a proper transfer function (of the electromagnetic wave) that is different than the Fourier transform. In this context, the most well studied case is the regime of Fresnel diffraction, where the transfer function is expressed in an integral form known as Fresnel integral . This regime occurs naturally at a range of distances away from the object plane which naturally includes also the Fraunhofer regime where the transfer function reduces to a simple Fourier transform. Going beyond the Fresnel regime is also possible. This means that the (magnitude squared of the) electromagnetic wave will be measured at some arbitrary plane away from the object. A more general case arisees by expressing the scalar transfer function of the light in a homogeneous medium, at any plane $z$ as:

Here $k=\omega/c$, with $\omega$ being the frequency of the light, $c$ being the speed of light in the medium, and $k_{x},k_{y}$ describe the transverse wavenumbers. The field at any arbitrary plane $z$, $E(x,y,z)$, is then given by inverse Fourier transforming the spectral function at that plane $F(k_{x},k_{y},z)$, which is related to the spectrum at the initial plane by:

With the transfer function 23, one can now formulate a new phase retrieval problem, where the measurements are conducted at some arbitrary plane $z$, giving $|E(x,y,z)|^{2}$, and the sought signal is $E(x,y,z=0)$. This approach can be extended to include polarization effects, where the transfer-function is vectorial, thereby describing the propagation through Maxwell’s equations with no approximation at all. The optical far-field - where the measurement corresponds to the Fourier magnitude of the image at the initial plane, (i.e., the measurement is proportional to $|F(k_{x},k_{y},z)|^{2}$) - is obtained for distances $z$ larger than some minimum distance $z_{0}$ that depends on the spectral extent of $F(k_{x},k_{y},z=0)$, and only within a region close enough to the $z$ axis in the measurement plane.

It is interesting to compare these more general phase retrieval problems to the generic problem of recovering a signal from its Fourier magnitude. In terms of algorithmics, the generic problem is much simpler and was extensively studied throughout the years, whereas the general case is considerably more complex and was studied only sporadically. However, in terms of optics, the measurements in the general case can provide more information. Namely, in the general case, measurements of $|E(x,y,z)|^{2}$ can be taken at multiple planes (multiple values of $z$), and each measurement adds more information on the signal. In contrast, for the generic problem, once the measurements are taken in the optical far-field, taking more measurements at further away distances does not add additional information because all the far field measurements correspond to the Fourier magnitude (to within some known scaling of coordinates in the measurement planes). As such, performing phase retrieval of optical images in the most general (non-Fourier) case can be beneficial, as it leads to multiple measurements, possibly relaxing the conditions on oversampling and/or the advance knowledge on the support in the image plane.

Historically, these ideas on non-Fourier measurements are known to the optics community since the early days of optical phase retrieval . They are currently being used in the context of improving the convergence of phase retrieval by taking non-Fourier measurements at several planes . Alternatively, one can take measurements at several different optical frequencies $\omega$, which would be expressed as different values of $k=\omega/c$ in the general transfer function given above. In this mutli-frequency context, it is important that the frequencies are well separated, each having a narrow bandwidth, to conform the high degree of coherence required for CDI. These ideas are now being pursued by several groups . Interestingly, the multi-frequency idea also works in the continuous case of broad bandwidth pulses centered on a single frequency. In this case, the power spectrum of the pulse must be known in advance . In a similar vein, recent work has demonstrated scanning CDI, where the beam is scanned through overlapping regions on the sample to allow imaging of extended objects, a method known as ptychography .

More sophisticated physical settings also exist, where the medium within which the waves are propagating is not homogeneous in space. Famous examples are photonic crystals, wherein the refractive index varies periodically in space, in a known fashion, in one, two or three dimensions. Obviously, in such settings the transfer function for electromagnetic waves is fundamentally different from the transfer function in free space. The phase retrieval problem in such systems, albeit less commonly known, is no less important. For example, photonic crystal fibers can in principle be used for imaging in endoscopy. The measurements in such systems correspond to the magnitude squared of the field at the measurement plane, which would be very different than the Fourier magnitude of the sought image. Still, once the transfer function is known, complicated as it may be, the phase retrieval problem is well defined and can be solved with some modifications to the algorithms described above. See, for example, pioneering work on phase retrieval in a photonic crystal fiber , and very recently on sparsity-based phase retrieval and super-resolution in optical waveguide arrays . In addition to these, the concept of CDI has also been extended to other schemes, such as Bragg CDI, suitable to periodic images, to reconstruct the structure and strain of nano-crystals .

VI-B Phase Retrieval Combining Holographic Methods

As explained in the introduction, optical settings always suffer from the inability of photodetectors to directly measure the phase of an electromagnetic wave at frequencies of THz (terahertz) and higher. Partial solution for this problem is provided through holography, invented by Denis Gabor in 1948 and awarded the Nobel Prize in Physics in 1971. Holography involves interfering an electromagnetic field carrying some image, $E_{image}$, with another electromagnetic field of the same frequency and a known structure, denoted as $E_{ref}$. Typically, the so-called reference wave, $E_{ref}$, has a very simply structure, for example, approximately a plane wave (wave of constant amplitude and phase). The detection system records $|E_{image}+E_{ref}|^{2}$. Originally, such holographic recording was done on a photographic plate which is made from a photosensitive material whose transmission becomes proportional to the recorded pattern $|E_{image}+E_{ref}|^{2}$. This photographic plate is called a hologram, wherein the information contained in the image wave $E_{image}$ is embedded in transmission function of the hologram. To see the recording, the wave of the known pattern, $E_{ref}$, is generated (which is possible because its structure is simple and fully known) and made to illuminate the hologram. The magnitude of the wave transmitted through the illuminated hologram is therefore proportional to $|E_{image}+E_{ref}|^{2}\cdot E_{ref}$. One of the terms is therefore $|E_{ref}|^{2}\cdot E_{image}$. Since $|E_{ref}|^{2}$ carries virtually no information (i.e., it is just a constant), this transmitted wave reconstructs the image times that constant. This is the principle of operation of holography. Over the years, it has been shown that it is almost always beneficial to record not the actual image but its Fourier spectrum, hence the reconstructed information is the Fourier transform of the image, and the image itself is recovered either in the far-field (as explained in the introduction) or at the focal plane of a lens. This process is termed Fourier holography .

In the context of phase retrieval, holography is used for the purpose of adding information in the measurement scheme. Because in most cases the measurements used are Fourier magnitudes, which physically implies far-field measurements, the natural inclusion of holographic methods is through Fourier holography. For example, adding a tiny hole (a delta function) at a predetermined position in the sample (close to where the sought image resides) creates an additional wave in the far field, with a tilted phase (arising from the displacement between the hole and the sought image). The far field intensity therefore now corresponds to the absolute value squared of the sum of the Fourier spectrum of the sought image and the (known) wave. As such, it introduces additional prior knowledge, which can be used for increased resolution of the algorithmic recovery or for relaxing the constraints on the prior knowledge on the support. These ideas have been exploited successfully using X-rays and electrons by several groups .

VI-C Challenges, Bottlenecks and Vision

The current challenges can be briefly defined in a single sentence: higher resolution, ability to recover more complex objects, improved robustness to noise, and real-time operation. The very reason phase retrieval in optical imaging has become so important is owing to the vision to be able one day to image complex biological molecules directly, track their structural evolution as it evolves in time, and even view the dynamics of the electronic wave functions bonding atoms together. The reasoning is obvious: to understand biology at the molecular level, to decipher the secrets of how their atomic constituents bond together and how they interact with other molecules. The current state of the art is far from those goals: imaging resolution is not yet at the atomic (sub-nanometer) level, and - at nanometric resolution - imaging cannot handle objects that are bounded by a support that is extremely large compared to the resolution. In terms of being able to perform real-time experiments, state of the art measurements have demonstrated extremely short optical pulses: tens of attoseconds ($10^{-18}$ seconds - on the order of the passage of a photon through a distance comparable to the size of an atom). Pioneering experiments have even started to probe the dynamics of electrons in molecules and tunneling processes on these time scales. But, as of today, none of these ultrafast methods was applied to imaging of even a simple molecule, let alone complex biological structures.

Clearly, the underlying physics and engineering pose great challenges to meet these goals. Generating coherent radiation in the hard X-ray regime is still a major obstacle, often requiring very large enterprises such as the X-ray sources at the SLAC National Accelerator Laboratory. These facilities around the world are continuously improving their photon flux at shorter wavelengths, thereby constantly improving imaging resolution. The fundamental limits on the coherent X-ray flux possible with current methods (such as synchrotrons, XFELs , and the process of high harmonics generation ) are not even known. But the steady improvement does give hope for imaging at the atomic level in the near future. Taking the CDI techniques to the regime of attosecond science is an important challenge. These pulses are extremely short, hence their bandwidth is huge, so the coherent diffraction pattern is a superposition of their Fourier contents, which requires new algorithmic methods. As described above, these issues are currently being explored by several groups. But the problem is fundamentally more complicated, because the process of scattering of light by molecules at these short wavelengths and ultrashort timescales is not like passing light through a mask on which an image is imprinted. Rather, many issues related to light-matter interactions under these conditions are yet to be understood (e.g., tunneling ionization of atoms by laser pulses).

Finally, the long term vision must include imaging the dynamics within complex biological systems at the atomic level and in real time. But such systems are extremely complex to handle, in terms of details on many spatial and temporal scales simultaneously, in terms of the statistical nature and huge redundancy in the physical processes taking place within such complexes simultaneously, and even in terms of the quantum mechanics governing the dynamics at those scales. This is where the signal processing community can make a large impact, by devising new and original methods for recovering the information from experimental measurements. Clearly, the algorithms will have to be tailored to the specific physical settings.

References