1. Introduction

In a wide range of wave phenomena in nature one does not have direct access to the signal under study but only to its time-average, losing the phase information. This is the characteristic case in optics [1], X-ray fields [2–6], high harmonic generation [7–9] and electron diffraction [10], often in the context of imaging. Retrieval of this inaccessible phase can be achieved by interferometry and holography, where one source, having a trivial or a-priori known phase is interfered with the signal of interest. Interferometry converts the phase modulation of the signal into amplitude changes in the interference pattern via the cross term that mixes the amplitudes and phases of both signals. Whereas in a simple interferometer a single phase is measured, this idea can be extended to the measurement of multiple phases by spatial multiplexing, as in holographic imaging. For imaging applications, holographic phase measurements can be performed either in the spatial domain or in the Fourier conjugate frequency domain [11]. Both of these have been extensively utilized in various applications such as quantitative phase imaging [7, 12] and multidimensional spectroscopy [13].

Yet, there exist many important physical scenarios where interferometry and holography are not feasible, particularly due to the difficulty in generating a stable, well-characterized reference beam. In such cases, various indirect pathways to retrieve the phase have been developed, known as coherent diffraction imaging (CDI) [2]. In CDI, a coherent wave illuminates an object, and its far field diffraction pattern is measured. If the illuminated object has a finite support and the diffraction pattern is sampled at a rate sufficiently higher than the Nyquist frequency such that the number of independently measured intensities is more than the number of unknown variables of the object, the phase information is commonly uniquely encoded in the diffraction intensities [14, 15]. However, the uniqueness of the solution does not imply the existence of a simple algorithm to retrieve it. In general, finding this solution is known to be a computationally difficult problem [16]. In fact, some instances of the phase retrieval problem are known to be nondeterministic polynomial (NP)-complete [17]. Over the years several iterative phase retrieval algorithms have been derived [1, 18, 19]. Despite broad use and notable successes, these iterative algorithms are associated with computational issues such as the convergence, stagnation and robustness to noise.

Both approaches described above - holography and CDI - have been independently applied in numerous applications [3, 5, 11, 20, 21]. In this paper, through a combination of these two approaches, we propose and demonstrate double blind Fourier holography as a new tool for 2D phase retrieval. On the one hand, our proposed method takes advantage of the holographic approach by measuring two independent components of the signal of interest and their interference. On the other hand, due to the oversampling of the diffraction intensities, the compact support can be utilized to overcome the absence of a well-characterized reference signal which is essential in holography [11,12,22,23]. In our method, both signals can be viewed symmetrically as an object and reference since both are unknown and can be, for instance, two arbitrary objects. We thus call this setup double blind Fourier holography. In contrast to common iterative 2D phase retrieval schemes, which in essence try to solve a challenging set of quadratic equations, in our approach the compact support constraint is translated to a set of linear equations, enabling the use of standard linear algebra tools for its solution. Moreover, our scheme allows a dimensionality reduction. Rather than solving a 2D problem, we dissect it into multiple 1D problems. The 1D case, whose properties were investigated in our previous work [24, 25], offers significant computational savings and is presented in more details below. Thus, this approach provides a direct and robust solution to the phase retrieval problem, without the need for either an external known reference as in holography, or an iterative reconstruction procedure as in CDI.

Double blind Fourier holography can be realized in a broad range of experimental scenarios. The two components of the signal can be associated with various attributes that yield different, overlapping Fourier transforms, such as two polarization vectors, spin states, propagation directions or quantum paths. As a proof of concept, we experimentally demonstrate this new scheme in a lensless optical imaging setup. In this case, our method might be viewed as holography in which both the object and reference fields are unknown and are complementary regions of a single 2D object. The fields are retrieved computationally using the compact support information.

2. Lensless imaging: problem setup and related works

CDI is a potent tool that has found applications in physics, material sciences, nanoscience and biology [3, 5, 20, 21, 26–29], recently attracting a significant attention in the context of X-ray imaging. In CDI, an object is illuminated by a spatially coherent plane wave and its far-field oversampled diffraction pattern is measured. For simplicity, we consider the case where the diffraction pattern is proportional to the square of the 2D Fourier transform of the object.

In the basic approach to CDI, only a single diffraction pattern is acquired. In this case, most phase retrieval methods are based on a combination of the oversampling method [14] with iterative algorithms, such as the Gerchberg-Saxton algorithm [30] and variants thereof developed by Fienup and others [1, 7, 17, 18, 31]. These iterative algorithms have to deal with computational difficulties such as their potentially slow convergence, stagnation and robustness to noise.

In general, the computational difficulties associated with a single diffraction measurement can be alleviated by an increased experimental complexity. One example is Fourier holography [11, 32], where a known reference beam is used to directly deduce the missing Fourier-domain phase. This scheme can be difficult or even unfeasible to apply, as the reference arm has to be interferometrically stabilized relative to the observed signal. An alternative approach is ptychography, also known as scanning CDI [33]. Here, rather than measuring a single diffraction pattern, multiple parts of the object are repeatedly measured through a limiting aperture at different spatial positions. The significant spatial overlap between the different exposures yields additional consistency constraints that accelerate the convergence of the iterative reconstruction algorithm. Ptychography has been successful in various scenarios [10, 34] and has mitigated convergence and robustness issues. However, it is not applicable in certain scenarios, such as few small and separated objects (compared to the scanning aperture size). Moreover, unlike our method, the robustness of ptychography depends on the number of measured apertures (See appendix E for comparison to ptychography).

Other methods, which come with mathematical guarantees of convergence and noise-robustness [35, 36], are computationally intensive and require an even higher degree of experimental complexity, for example by illuminating the sample through multiple masks with known parameters and positions. Although promising, these have not yet been experimentally demonstrated.

3. The proposed method

Our approach is based on diffraction pattern measurements of two independent signals, and recovery of their interference. The interference measurement already reveals significant information about the relative phase between the signals. This is, however, insufficient for full reconstruction since multiplication of both measured signals by a common phase pattern will not alter the result of any of the measurements. We have previously shown that in the 1D case, the constraint of a (possibly unknown) compact support is sufficient to fully resolve the phases of both signals, as it constrains the common phase to a unique solution [24, 25]. Moreover, this solution can be obtained by solving a set of linear equations. This approach can be readily generalized to the 2D case. Yet, for a n × n image this would lead to a very large system of linear equations with n² unknowns, not amenable to direct inversion.

A natural pathway to address this difficulty is to dissect the 2D problem into a set of many independent 1D blind-holography problems and then properly stitch the 1D solutions to the 2D one. This provides a direct pathway for solving the 2D double blind holography problem. The reduced dimensionality is in general not possible in other phase retrieval schemes except for zeros localization methods as in [37], since in general there is no unique solution to the 1D phase problem. This reduction enables scalability of the reconstruction procedure to very large images. In the following, we first briefly describe the attributes of the 1D scheme (also termed VPR - vectorial phase retrieval), then show its generalization to two dimensions, and conclude with an experimental demonstration in the optical domain.

The basic principles for the 1D case are summarized below. For simplicity we focus on the noise-free measurement of a signal of length n with a known support of length n/2 (see appendices A and B for a detailed description and [25] for a more general case of noisy measurements of a signal whose compact support is unknown). Let us consider a two component vectorial signal E_1,2(x_i) of length n of which E_1,2(x_i) = 0 for i > n/2. Importantly, these two components need not be scalar components of a particular physical vector. ’Vectorial’ merely means that the two components are packed in a vector form. The Fourier Transform of the signal is represented as Ẽ_1,2(k) = |Ẽ_1,2(k)|e^iϕ_1,2(k) where k is the Fourier domain coordinate and ϕ_1,2(k) is the unknown phase in the Fourier domain. With |Ẽ_1,2(k)|² and their interference, ${\tilde{E}}_1(k){\tilde{E}}_2^{*}(k)$ as input, 1D VPR recovers the vector of 2n unknown phases X_1,2(k) = e^iϕ_1,2(k) based on the following set of linear equations. First, the interference measurements provide n linear equations,

One can envision numerous applications of the extension of 1D VPR to two dimensions. As an example, we focus on the case of lensless imaging. Such 2D phase retrieval problems have been a subject of intensive studies, and remain a challenge in practical lensless applications.

The realization of 2D VPR in a lensless imaging scenario is schematically illustrated in Fig. 1(a) in which an object is defined by a 2n_x × n_y aperture (Fig. 1(b)). As described above, conventional lens-less methods reconstruct the object by first measuring a single oversampled diffraction pattern and then applying an iterative algorithm. In 2D VPR, in contrast, we first split the image into two parts E_1,2(x, y), associated with the two smaller n_x × n_y apertures shown in Fig. 1(c). These define the two vectorial components of the image. We then follow the scheme illustrated in Fig. 2. The top panels in Fig. 2 show three 2D diffraction patterns, measured with the two small apertures (denoted as |Ẽ₁(k⃗)|² and |Ẽ₂(k⃗)|²) and the large aperture (denoted as |Ẽ₃(k⃗)|² = |Ẽ₁(k⃗) + Ẽ₂(k⃗)|²). The 2D VPR scheme first determines ${\tilde{E}}_1(\overrightarrow{k}){\tilde{E}}_2^{*}(\overrightarrow{k})$ from the three diffraction pattern measurements. The real part of this quantity is resolved by the following identity,

 

Fig. 1 Schematic layout of double blind Fourier holography in lensless imaging method. (a) Three diffraction patterns are measured from a sample by varying the illumination area. The varied illumination area is produced by a large (b) and two small (c) apertures. The apertures are positioned in proximity to the sample (either before the object or after it). The small apertures (b) are at two neighboring positions, their union covers the same area as the large aperture (c).

Download Full Size | PDF

 

Fig. 2 Image reconstruction based on double blind Fourier holography. The scheme is demonstrated using experimental data of a sample consisting of the letters WIS. The sample is illuminated by a Gaussian laser beam. Top panel: three diffraction patterns measured from the letter W (left), the letters IS (middle) and the entire WIS structure (right). Each diffraction pattern is sliced into columns (middle left) and rows (middle right). For clarity, only two columns and rows are shown. The rows and columns are then solved by the 1D-phase-retrieval algorithm to obtain spatial-spectral maps (bottom left and bottom right panels). The two reconstructed spatial-spectral maps are then stitched (see section 3) to form the reconstructed object (bottom middle panel). The relative intensity of the reconstructed letters WIS faithfully recovers the Gaussian illumination pattern used in the experiment.

Download Full Size | PDF

At this point, given |Ẽ₁(k⃗)|², |Ẽ₂(k⃗)|² and ${\tilde{E}}_1(\overrightarrow{k}){\tilde{E}}_2^{*}(\overrightarrow{k})$ we could apply our VPR scheme directly (by applying a 2D version of Eqs. (1) and (2)). This, however, leads to a very large system of linear equations with 𝒪(n²) unknowns. Instead, we now show how the 2D problem can be dissected into a set of many 1D problems. This can be done by applying the 1D VPR for each row and coloumn of the 2D signal separately, and then combining them together. To this end, we decompose |Ẽ₁(k⃗)|², |Ẽ₂(k⃗)|² and ${\tilde{E}}_1(\overrightarrow{k}){\tilde{E}}_2^{*}(\overrightarrow{k})$ into a set of 1D slices along the horizontal direction, as schematically shown in Fig. 2 (for graphical clarity we describe only 2 such slices, denoted a and b). The slices of the three diffraction patterns are then combined to form 1D vectorial signals. By using the 1D VPR scheme, we can directly retrieve the phase vector associated with each slice up to an arbitrary global phase e^iϕ(k_y) (for each 1-D slice), which depends on the slice index k_y.

We repeat the above process in the perpendicular direction and obtain the phase vectors associated with each vertical slice (denoted c and d in Fig. 2), up to an arbitrary global phase e^iψ(k_x), which depends on the slice index k_x. Finally, the 2D image is found by a consistent stitching of the vertical and horizontal solutions. The two solutions, which represent the same object, should be equal to each other up to a multiplicative factor e^{i(ψ(k_x)−ϕ(k_y))}. Their ratio is therefore a rank 1 matrix, whose left and right eigenvectors are e^iψ(k_x) and e^−iϕ(k_y). These can be easily found numerically by a singular value decomposition (SVD) and subtracted from the reconstruction to give the correct 2D phase, which is consistent with the compact support in both directions. For more details on the algorithm see appendices A–D.

4. Experimental demonstration

To illustrate the 2D double blind holography scheme, we built an experimental system based on a He-Ne laser with a wavelength of 632.8 nm. The sample is a transmission mask milled onto a thin metallic plate, shown in Fig. 3(e). Three diffraction patterns were acquired from the sample by using two small apertures and a larger aperture. The illumination of the two small apertures was achieved by constructing a movable slit to define two regions within the sample. The diffraction patterns were measured by a charge-coupled device (CCD) at a resolution providing sufficient oversampling for direct phase determination. Figures 3(a)–(c) show the three diffraction patterns, acquired with two small apertures and a large one, respectively. By using the 2D phase retrieval scheme discussed above, we directly determined the phase information of the diffraction intensity and reconstructed the sample structure (Fig. 3(d)), which is in good agreement with the direct image of the sample (Fig. 3(e)). In our direct phase determination, the size and position of the compact support were unknown and computationally identified by our scheme (see appendix D for more details). In the experiment the object contained a zero valued region which we used to define the boundary between the two different components. If no such region exists, the experimental difficulty increases due to the necessity to prevent overlap between the two components. However, this is the only major alignment challenge, since as long as there is no overlap, our method is oblivious to the exact shape and position of the compact support.

 

Fig. 3 Experimental demonstration of double blind Fourier holography in lensless imaging. The sample (e) comprises of a transmission mask of a tree, milled onto a thin metallic plate (the diameter of the circular sample is 3.5mm) The vertical black stripe in the middle (which is a part of the image) defines the two apertures. Panels a–c: Three diffraction patterns measured from the two small and one large aperture, respectively. The two small apertures were realized by constructing a movable slit to define the left and right parts of the sample. d, 2D image reconstruction from the three diffraction patterns using 2D VPR. The reconstructed image is in good agreement with the original sample structure.

Download Full Size | PDF

5. Conclusions

In conclusion, in this manuscript we propose and experimentally demonstrate a conceptually new scheme for direct 2D phase determination, based on only three diffraction pattern measurements. Our method reduces the 2D phase problem to a set of linear equations and is hence scalable, efficient, and robust to noise. While the proof-of-principle experiment was conducted with an optical laser, this method can be readily adapted to the reconstruction of complex objects using other sources such as high harmonic generation and synchrotron radiation. Although we focused on the lensless imaging application in this manuscript, due to the ubiquitous nature of the phase problem the proposed double blind Fourier holography scheme might have applications in other fields as well.

Appendix A:general setup

The inputs of our proposed method are three measurements of diffraction patterns of an object with a finite support: the first two, |Ẽ₁(k⃗)|² and |Ẽ₂(k⃗)|², are both taken with a small aperture at neighboring positions, aligned in the x axis. The third pattern, |Ẽ₃(k⃗)|² = |Ẽ₁(k⃗) + Ẽ₂(k⃗)|², is taken with a large aperture of double size in the x direction (see Fig. 1(b) and (c)). Under the Fraunhofer approximation, each of the diffraction patterns Ẽ_1,2,3(k⃗) is given by the 2D discrete Fourier transform (DFT) of the corresponding E_1,2,3(x⃗). It is useful to recall that the 2D DFT ℱ_2D can be represented as a sequence of two 1D DFT’s ℱ_x and ℱ_y in the x and y directions, respectively:

Let the small aperture be of size n_x × n_y pixels and the larger aperture be of size 2n_x × n_y. For simplicity, we consider sampling of the diffraction pattern at exactly the Nyquist rate of E₃ (the advantages of denser sampling will be discussed later), so that the measured intensities |Ẽ_1,2,3(k⃗)|² are images of size 4n_x × 2n_y. These measured patterns are related to the original signals E_1,2,3(x⃗) through a zero padded 2D DFT. Since we measure absolute values, the exact location of the padding can be chosen arbitrarily, as it only introduces a multiplicative phase factor.

The apertures implicitly introduce a finite support on E_1,2,3(x⃗), of size at most n_x × n_y for E₁, E₂ and 2n_x × n_y for E₃. Given the arbitrary choice of location for the zero-padding, without loss of generality, we assume that in the spatial domain the three signal’s supports are given by:

(4)$$\begin{array}{lll}{S}_1\hfill & =\hfill & \left\{x\le n_x\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}y\le n_y\right\}\hfill \\ S_2\hfill & =\hfill & \left\{n_x<x\le 2n_x\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}y\le n_y\right\}\hfill \\ S_3\hfill & =\hfill & \left\{x\le 2n_x\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\le n_y\right\}.\hfill \end{array}$$

Next, note that even though E_j(x⃗) is supported on the subset of indices S_j and is zero elsewhere, its 2D Fourier transform Ẽ_j(k⃗) is in general non-zero at all 4n_x × 2n_y pixels. The key principle that allows dissection of the 2D VPR problem into a sequence of 1D problems is based on the following observation: the intermediate 1D Fourier transforms, ℱ_xE_j and ℱ_yE_j, do have a finite support (see Fig. 4). In particular, E_1,2(x, k_y) are zero for x > n_x and E₃(x, k_y) is zero for x > 2n_x. Similarly, E_1,2,3(k_x, y) are zero for y > n_y. This property, illustrated in Fig. 4, will be useful in the developments below.

 

Fig. 4 Fourier transform changes the support in the direction it acts, but maintain the support in the perpendicular direction.

Download Full Size | PDF

Appendix B: reconstructing ${\tilde{E}}_1{\tilde{E}}_2^{*}$ from |Ẽ₁|, |Ẽ₂|, |Ẽ₃|

Our vectorial phase retrieval scheme requires, as input, both |Ẽ₁|², |Ẽ₂|² and the interference ${\tilde{E}}_1(\overrightarrow{k}){\tilde{E}}_2^{*}(\overrightarrow{k})$. The first two are directly measured. As for the third one, its real part can be easily recovered from the diffraction measurements, since, as stated in Eq. (3):

$$\mathit{Re}\left[{\tilde{E}}_1{\tilde{E}}_2^{*}\right]=\frac{1}{2}\left({\left|{\tilde{E}}_3\right|}^2-{\left|{\tilde{E}}_1\right|}^2-{\left|{\tilde{E}}_2\right|}^2\right).$$

In our previous works [24, 25] we assumed a fourth measurement |Ẽ₄| = |Ẽ₁ + iẼ₂|, which allowed recovery of the imaginary part by a formula similar to Eq. (3).

Here, we show how to recover the imaginary part without such an additional measurement, but at the expense of sampling the Fourier transform at a higher rate. Indeed, this factor 2 sampling rate is implicit in our description, since we assumed measurements of 4n_x Fourier coefficients in the x direction (instead of 2n_x, which is the Nyquist rate for E₁, E₂).

To recover the required imaginary part, note that by the convolution theorem for 1D Fourier transforms,

$${\tilde{E}}_1(\overrightarrow{k}){\tilde{E}}_2^{*}(\overrightarrow{k})={ℱ}_x\left[E_1(x,k_y){*}_x{E}_2(4n_x-x,k_y)\right],$$

where *_x denotes convolution in the x direction. As E₁(x, k_y) vanishes for all x > n_x and E₂(4n_x − x, k_y) vanishes for all x < 3n_x, their convolution E₁(x, k_y) *_x E₂(4n_x − x, k_y) has a support of size at most 2n_x − 1. Moreover, its support is not arbitrary but rather given by the indices x = 2,...,2n_x.

Let us denote by $R(x,k_y)={ℱ}_x^{-1}\left[\mathit{Re}\left\{{\tilde{E}}_1(\overrightarrow{k})E_2^{*}(\overrightarrow{k})\right\}\right]$ and $I(x,k_y)={ℱ}_x^{-1}\left[\mathit{Im}\left\{{\tilde{E}}_1(\overrightarrow{k})E_2^{*}(\overrightarrow{k})\right\}\right]$. For any fixed k_y, both R(x, k_y) and I(x, k_y) are vectors of length 4n_x. From Eq. (5) it follows that R(x, k_y) is known and hence our goal is to recover I(x, k_y). By the linearity of the DFT and the fact that E₁(x, k_y) *_x E₂(−x, k_y) is non-zero only at the indices {2,...,2n_x}, we have

(5)$$R(x,k_y)+\mathbf{i}I(x,k_y)=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{for}\hspace{0.17em}\hspace{0.17em}x=1\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\text{for}\hspace{0.17em}\hspace{0.17em}x=2n_x+1,\dots ,4n_x.$$

That is, I(x, k_y) = iR(x, k_y) at these indices. Finally, since I(x, k_y) is the inverse DFT of a real valued function,

$$I(x,k_y)=I^{*}(4n_x-x,k_y).$$

Together with Eq. (5) this allows complete reconstruction of I(x, k_y) from R(x, k_y).

Remark: The above scheme is intimately related to the Kramers-Kronig relations [38], whereby the imaginary part of an analytic signal is recovered from knowledge of its real part. In our setup the support of E₁(x, k_y) *_x E₂(x, k_y) plays a role analogous to that of temporal causality in the Kramers-Kronig case.

Appendix C: the 2D VPR reconstruction algorithm

As outlined in the main text, our proposed algorithm consists of three steps:

1. Reconstruct ${\tilde{E}}_1(\overrightarrow{k}){\tilde{E}}_2^{*}(\overrightarrow{k})$ from |Ẽ₁|, |Ẽ₂|, |Ẽ₃|;
2. Reconstruct the phases independently for each row and column of the 2D signal. At this point we remain with an overall phase for each row and column which we denote by e^iϕ(k_y) and e^iψ(k_x).
3. Reconstruct, up to a single global phase, the unknown phases e^iϕ(k_y) and e^iψ(k_x).

Step i) was described in the previous section. For a more detailed description of the two remaining steps, we now introduce some notations. First, recall our notation for the unknown phases, Ẽ_1,2(k⃗) = |Ẽ_1,2(k⃗)|X_1,2(k⃗), where X_1,2(k) = e^iϕ_1,2(k) are the unknown phases to be retrieved. As discussed in the main text, for each fixed k_y = m, we obtain a 1-D vectorial phase retrieval problem. For each fixed k_y = m we combine the two unknown vectors of phases X₁(k_x, m) and X₂(k_x, m) into a single unknown vector of length 8n_x,

$${\text{x}}_m=\left[X_1(k_x,m),X_2(k_x,m)\right].$$

Following our previous works [24, 25], the vector x_m can be recovered via the non-trivial solution to the following set of homogeneous linear equations

$$A_x{\text{x}}_m=0,$$

where the matrix A_m is of size 8n_x × 8n_x, and its entries are given by

(6)$${\left(A_m\right)}_{i,j}=\{\begin{array}{ll}\left|{\tilde{E}}_1(j,m)\right|e^{\mathbf{i}\frac{j(i+2n_x)2\pi }{4n_x}}\hfill & i\le 2n_x,j\le 2n_x\hfill \\ \left|{\tilde{E}}_2(j-2n_x,m)\right|e^{\mathbf{i}\frac{(j-2n_x)i2\pi }{4n_x}}\hfill & 4n_x\ge i>2n_x,4n_x\ge j>2n_x\hfill \\ {\tilde{E}}_1(j-4n_x,m){\tilde{E}}_2^{*}(j-4n_x,m)\hfill & j>4n_x,j-4n_x=i\hfill \\ -\left|{\tilde{E}}_1(j-4n_x,m)\right|\left|{\tilde{E}}_2(j-4n_x,m)\right|\hfill & j>4n_x,j=i\hfill \\ 0\hfill & \mathit{otherwise}\hfill \end{array}$$

As proven in [25], for spectrally independent signals in the ideal case of noise-free measurements, the matrix A_m has a singular value 0 with multiplicity one, whose corresponding singular vector (upon proper normalization) is a phase vector yielding the exact solution to the phase retrieval problem.

In the presence of noise, we look for the column of A_m with maximal norm, and set the corresponding entry of x_m to one. We denote the vector x_m without this index by x_m−, the corresponding column of A_m by $a_m^1$ and the matrix A_m without this column by (A_m)₋. For each k_y = m, we then compute the minimizer x̂_m of ${\Vert (A_m)-{\text{x}}_{m-}-a_m^1\Vert }^2$. This recovers x_m, up to an unknown multiplicative phase e^iϕ(k_y). If the data is over-sampled at a rate higher than the Nyquist rate, then there are more equations than unknowns, which improves the robustness to noise.

Next, we repeat the above procedure in the perpendicular direction, solving for k_x = m the unknown phases y_m = [X₁(m, k_y), X₂(m, k_y)] up to an unknown multiplicative phase e^iψ(x).

At this point, since the 1-D VPR was applied independently to each row and column, we are left with an overall unknown phase for each one, denoted by e^iϕ(k_y) and e^iψ(k_x). The last step in reconstructing the 2D object E_1,2(x, y), is to recover the those phases. To this end, we construct the matrix $z(k_x,k_y)\frac{\text{x}(k_x,k_y)}{\text{y}(k_x,k_y)}$ where x(k_x, k_y) and y(k_x, k_y) are the phases reconstructed in the previous step (by applying 1-D VPR) for the rows and the columns respectively, and the ratio is defined element-wise. This is implemented by multiplying x(k_x, k_y) by the complex conjugate of y(k_x, k_y) instead of using explicit division. Since in the noise free case both x and y contain the exact phase of Ẽ_1,2(k⃗) up to a multiplicative phase factor e^iψ(k_x) or e^iϕ(k_y), their ratio is of the form:

$$z(k_x,k_y)=\frac{\text{x}(k_x,k_y)}{\text{y}(k_x,k_y)}=e^{\mathbf{i}\left(\varphi (k_y)-\psi (k_x)\right)}.$$

Namely, z(k_x, k_y) is a rank 1 matrix, whose left and right singular vectors are e^iψ(x) and e^−iϕ(y). Via an SVD of this matrix we thus recover the unknown phase factors e^iψ(x) and e^−iϕ(y). Subtracting these from x(k_x, k_y) and y(k_x, k_y) respectively gives the solution of the phase problem.

Table 1 summarizes the algorithm for solving the 2D phase problem with the proposed experimental setup.

Table 1:. Algorithm for solving the vectorial lensless phase problem

View Table

Appendix D: unknown compact support

The algorithm described above assumes a known compact support, which is not always the case. To deal with this issue, the 1-D vectorial phase retrieval algorithm in fact contains a model selection step which estimates the compact support, as described in [25]. The basic idea behind this step is that while the system of linear equations we solve always has some solution (in least squares sense), for incorrect assumed compact support, the matrix A_m either has a minimal singular value strictly larger than zero, or has multiplicity at the zero singular value. In both cases, the resulting solution is not a phase vector and upon normalization into a phase vector it has a larger residual energy outside the assumed support. Thus, by scanning the possible supports length and checking the residual outside the assumed support we can estimate the correct support length. Note that our scan for the support size is not a 2D scan, as we can look for the support length separately for the rows and for the columns.

Figure 5 shows the results of the compact support model selection scans for the data shown in figure 3 of the main text. The x axis represent the assumed support length and the y axis represent the normalized energy of the VPR solution outside the assumed support, as described in [25]. The estimated compact support is the x value for which the energy outside the compact support is minimal

 

Fig. 5 Estimating the compact support: The right panel shows a scan over the support of the horizontal axis. The left panel shows a scan over the assumed support of the small aperture in the vertical axis.

Download Full Size | PDF

Appendix E: comparison to Ptychography

Ptychography is a lensless imaging method in which the image is reconstructed from a series of diffraction patterns, each taken with the same aperture at a different location. The spatially overlapping regions between the shifted aperture locations are used as additional constraints in the phase retrieval procedure. The image is reconstructed using an iterative reconstruction method, known as PIE (Ptychographic Iterative Engine). The experimental similarity between PIE and double blind Fourier holography (DBFH) is obvious: both methods require illumination of the sample with multiple different (or shifted) apertures. However, the working principle is quite different as are several computational and experimental aspects. In the following we list some of the similarities and differences between the methods, followed by numerical reconstructions for the noise-free and noisy cases for comparison.

1. Principle of operation: PIE is based on Fienup-type iterative algorithms utilizing compact support and additional constraints from the multiple aperture positions (and possibly other object dependent constraints). In contrast, DBFH is based on a holographic scheme without a known reference field, and compensates for the absence of a known reference by using the compact support constraint. As a consequence, while PIE iteratively solves a nonlinear problem, DBFH uses the interference term to reduce it to a linear one, solved by linear algebra tools.
2. Convergence: PIE has no mathematical guarantees of convergence to the exact solution even for the noise-free measurements (even though in practice it does converge for many cases). In contrast, DBFH is based on solving a set of linear equations, which is guaranteed to yield the correct solution if it is unique, which is the general case [39].
3. Apertures positioning: PIE requires one aperture while DBFH requires two apertures of different dimensions. While PIE is sensitive to the exact positioning of the apertures [40] DBFH only requires that there is no overlap between the small apertures. Moreover, if there is a zero valued region within the objects (or simply two separated objects), as in this paper, this potential difficulty of DBFH is also alleviated.
4. Number of apertures: PIE depends on the number and overlap of the apertures (to get sufficient additional constraints). In fact, for the object below, PIE did not converge to the correct solution using only 3 aperture positions. In contrast, DBFH requires exactly 3 apertures.
5. Small objects: PIE cannot be used for objects which are significantly smaller than the mask. In contrast, for DBFH it is possible to image two separate small objects, measuring the diffraction from each of them separately and from both of them together.
6. Numerical experiment: In practice, as numerically demonstrated below, PIE converges very quickly to the correct solution, almost irrespectively of the choice of the initial guess, for some objects. Yet, for other objects, it stagnates at local minima and almost never reconstructs the object correctly. Exact conditions for convergence of ptychography are not known. For example, as can be seen below, for the object reconstructed in this paper (see Fig. 3) PIE does not converge to the correct solution even in the noise-free case (for 5 aperture positions with 50% overlap). In contrast, for noise-free measurements DBFH always reconstructs the object if the problem has a unique solution.

To compare noiseless reconstructions, we considered two different objects: the Weizmann tree mask on which we performed our experiment (object 1) and the sample (object 2) proposed in http://www.psi.ch/sls/csaxs/PublicationsEN/thesis_dierolf.pdf, which is also the source for the PIE matlab code we used for the simulations. In both cases, PIE does not converge to the correct solution with 3 aperture positions used, but converges to a reasonable solution for 5 positions with 50% overlap. For object 2, PIE always converged to nearly the exact solution after about 150 iterations (Fig. 6). In contrast, it never converged to the exact solution for object 1, though the reconstruction is similar to the original object (Fig. 7). DBFH reconstructed both objects exactly (with an error comparable to the machine error). A detailed study of the reason why PIE does not converge for some objects is beyond the scope of this work.

 

Fig. 6 PIE reconstruction of object 2 with 5 aperture positions and 50% overlap: Left -original object. Right - its reconstruction after 500 iterations. The reconstruction is nearly perfect.

Download Full Size | PDF

 

Fig. 7 PIE reconstruction of object 1 with 5 aperutre positions and 50% overlap: Left -original object. Right - its reconstruction after 1000 iterations. Bottom - Error vs. iteration number. As can clearly be seen, the algorithm stagnates at a non-zero error, and the reconstruction is not perfect even for the noise-free case.

Download Full Size | PDF

In order to compare robustness to noise, we performed a Monte-Carlo simulation, in which we used noisy diffraction patterns to reconstruct the objects, and plotted the median of the reconstruction mean squared error (MSE) as a function of the noise level. The median is plotted rather than the average, since there exist rare events in which the PIE doesn’t converge to a reasonable reconstruction. These bias the average but not the median. 100 realizations of an uncorrelated zero mean Gaussian noise were used. For PIE we used 5 aperture positions with 50% overlap. We performed the numerical simulations for both object 1 and object 2, as can be seen in Fig. 8. Since PIE uses 5 positions and DBFH 3 aperture positions, the overall signal to noise ratio in PIE is higher. This explains the differences in the reconstruction error at low noise level of object 2.

 

Fig. 8 Reconstruction error vs. noise level: Left - object 1, Right - Object 2.

Download Full Size | PDF

Acknowledgments

D.O. and N.D. acknowledge support from the Israeli Centers of Research Excellence program and the Crown photonics center. N.D. acknowledges support by the Israeli Science Foundation and the Minerva Foundation. N.D. acknowledges support by the European Research Council starting investigator grant MIDAS. D.O. acknowledges support by the European Research Council starting investigator grant SINSLIM 258221. B.N. acknowledges support by the Israeli Science Foundation. J.M. acknowledges the support by the DARPA PULSE program through a grant from AMRDEC.

References and links

1. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef]   [PubMed]  

2. J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of x-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature 400, 342–344 (1999). [CrossRef]  

3. H. Jianga, C. Songb, C. C. Chena, R. Xua, K. S. Rainesa, B. P. Fahimiana, C. H. Luc, T.-K. Leec, A. Nakashimad, J. Uranod, T. Ishikawab, F. Tamanoi, and J. Miao, “Quantitative 3d imaging of whole, unstained cells by using x-ray diffraction microscopy,” Proc. Natl. Acad. Sci. USA 107, 11234–11239 (2010). [CrossRef]  

4. H. N. Chapman, P. Fromme, A. Barty, T. A. White, R. A. Kirian, A. Aquila, M. S. Hunter, J. Schulz, D. P. DePonte, U. Weierstall, R. B. Doak, F. R. N. C. Maia, A. V. Martin, I. Schlichting, L. Lomb, N. C. ans, R. L. Shoeman, S. W. Epp, R. Hartmann, D. Rolles, A. Rudenko, L. Foucar, N. Kimmel, G. Weidenspointner, P. Holl, M. Liang, M. Barthelmess, C. Caleman, S. Boutet, M. J. Bogan, J. Krzywinski, C. Bostedt, S. Bajt, L. Gumprecht, B. Rudek, B. Erk, C. Schmidt, A. Hmke, C. Reich, D. Pietschner, L. Strder, G. Hauser, H. Gorke, J. Ullrich, S. Herrmann, G. Schaller, F. Schopper, H. Soltau, K. U. Khnel, M. Messerschmidt, J. D. Bozek, S. P. Hau-Riege, M. Frank, C. Y. Hampton, R. G. Sierra, D. Starodub, G. J. Williams, J. Hajdu, N. Timneanu, M. M. Seibert, J. Andreasson, A. Rocker, O. Jnsson, M. Svenda, S. Stern, K. Nass, R. Andritschke, C. D. Schrter, F. Krasniqi, M. Bott, K. E. Schmidt, X. Wang, I. Grotjohann, J. M. Holton, T. R. Barends, R. Neutze, S. Marchesini, R. Fromme, S. Schorb, D. Rupp, M. Adolph, T. Gorkhover, I. Andersson, H. Hirsemann, G. Potdevin, H. Graafsma, B. Nilsson, and J. C. Spence, “Femtosecond x-ray protein nanocrystallography,” Nature 470, 73–77 (2011). [CrossRef]   [PubMed]  

5. M. M. Seibert, T. Ekeberg, F. R. N. C. Maia, M. Svenda, J. Andreasson, O. Jonsson, D. Odic, B. Iwan, A. Rocker, D. Westphal, M. Hantke, D. P. DePonte, A. Barty, J. Schulz, L. Gumprecht, N. Coppola, A. Aquila, M. Liang, T. A. White, A. Martin, C. Caleman, S. Stern, C. Aberge, V. Seltzer, J.-M. Claverie, C. Bostedt, J. D. Bozek, S. Boutet, A. A. Miahnahr, M. Messerschmidt, J. Krzywinski, G. Williams, K. O. Hodgson, M. J. Bogan, C. Y. Hampton, R. G. Sierra, D. Starodub, I. Andersson, S. Bajt, M. Barthelmess, J. C. H. Spence, P. Fromme, U. Weierstall, R. Kirian, M. Hunter, R. B. Doak, S. Marchesini, S. P. Hau-Riege, M. Fran, R. L. Shoeman, L. Lomb, S. W. Epp, R. Hartmann, D. Rolles, A. Rudenko, C. Schmidt, L. Foucar, N. Kimmel, P. Holl, B. Rudek, B. Erk, A. Hmke, C. Reich, D. Pietschner, G. Weidenspointner, L. Strder, G. Hauser, H. Gorke, J. Ullrich, I. Schlichting, S. Herrmann, G. Schaller, F. Schopper, H. Soltau, K.-U. Khnel, R. Andritschke, C.-D. Schrter, F. Krasniqi, M. Bott, S. Schorb, D. Rupp, M. Adolph, T. Gorkhover, H. Hirsemann, G. Potdevin, H. Graafsma, B. Nilsson, H. N. Chapman, and J. Hajdu, “Single mimivirus particles intercepted and imaged with an x-ray laser,” Nature 470, 78–81 (2011). [CrossRef]   [PubMed]  

6. H. N. Chapman and K. A. Nugent, “Coherent lensless x-ray imaging,” Nat. Photon. 4, 833–839 (2010). [CrossRef]  

7. R. L. Sandberg, C. Song, P. W. Wachulak, D. A. Raymondson, A. Paul, B. Amirbekian, E. Lee, A. E. Sakdinawat, L.-O. Chan, M. C. Marconi, C. S. Menoni, and M. M. Murnane, “High numerical aperture tabletop soft x-ray diffraction microscopy with 70-nm resolution,” Proc. Natl. Acad. Sci. USA 105, 24–27 (2008). [CrossRef]  

8. A. Ravasio, D. Gauthier, F. R. N. C. Maia, M. Billon, J. P. Caumes, D. Garzella, M. Géléoc, O. Gobert, J. F. Hergott, A. M. Pena, H. Perez, B. Carré, E. Bourhis, J. Gierak, A. Madouri, D. Mailly, B. Schiedt, M. Fajardo, J. Gautier, P. Zeitoun, P. H. Bucksbaum, J. Hajdu, and H. Merdji, “Single-shot diffractive imaging with a table-top femtosecond soft x-ray laser-harmonics source,” Phys. Rev. Lett. 103, 028104 (2009). [CrossRef]   [PubMed]  

9. R. L. Sandberg, A. Paul, D. A. Raymondson, S. Hädrich, D. M. Gaudiosi, J. Holtsnider, R. I. Tobey, O. Cohen, M. M. Murnane, H. C. Kapteyn, C. Song, J. Miao, Y. Liu, and F. Salmassi, “Lensless diffractive imaging using tabletop coherent high-harmonic soft-x-ray beams,” Phys. Rev. Lett. 99, 098103 (2007). [CrossRef]   [PubMed]  

10. M. Humphry, B. Kraus, A. Hurst, A. Maiden, and J. Rodenburg, “Ptychographic electron microscopy using high-angle dark-field scattering for sub-nanometre resolution imaging,” Nat. Commun. 3, 730 (2012). [CrossRef]   [PubMed]  

11. S. Eisebitt, J. Lüning, W. Schlotter, M. Lörgen, O. Hellwig, W. Eberhardt, and J. Stöhr, “Lensless imaging of magnetic nanostructures by x-ray spectro-holography,” Nature 432, 885–888 (2004). [CrossRef]   [PubMed]  

12. M. Mir, B. Bhaduri, R. Wang, R. Zhu, and G. Popescu, Quantitative Phase Imaging (Progress in Optics: Elsevier Science Limited, 2012).

13. S. Mukamel, “Multidimensional femtosecond correlation spectroscopies of electronic and vibrational excitations,” Annu. Rev. Phys. Chem. 51, 691–729 (2000). [CrossRef]   [PubMed]  

14. J. Miao, D. Sayre, and H. Chapman, “Phase retrieval from the magnitude of the fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A 15, 1662–1669 (1998). [CrossRef]  

15. J. Miao, T. Ishikawa, E. H. Anderson, and K. O. Hodgson, “Phase retrieval of diffraction patterns from noncrystalline samples using the oversampling method,” Phys. Rev. B 67, 174104 (2003). [CrossRef]  

16. S. Marchesini, “A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007). [CrossRef]  

17. H. Sahinoglou and S. D. Cabrera, “On phase retrieval of finite-length sequences using the initial time sample,” IEEE Trans. CIRC Sys. pp. 954–958 (1991). [CrossRef]  

18. S. Marchesini, H. He, H. N. Chapman, S. P. Hau-Riege, A. Noy, M. R. Howells, U. Weierstall, and J. C. H. Spence, “X-ray image reconstruction from a diffraction pattern alone,” Phys. Rev. B 68, 140101 (2003). [CrossRef]  

19. J. A. Rodriguez, R. Xu, C.-C. Chen, Y. Zou, and J. Miao, “Oversampling smoothness: an effective algorithm for phase retrieval of noisy diffraction intensities,” J. Appl. Crystallogr. 46, 312–318 (2013). [CrossRef]   [PubMed]  

20. X. Huang, J. Nelson, J. Kirz, E. Lima, S. Marchesini, H. Miao, A. M. Neiman, D. Shapiro, J. Steinbrener, A. Stewart, J. J. Turner, and C. Jacobsen, “Soft x-ray diffraction microscopy of a frozen hydrated yeast cell,” Phys. Rev. Lett. 103, 198101 (2009). [CrossRef]  

21. Y. Nishino, Y. Takahashi, N. Imamoto, T. Ishikawa, and K. Maeshima, “Three-dimensional visualization of a human chromosome using coherent x-ray diffraction,” Phys. Rev. Lett. 102, 018101 (2009). [CrossRef]   [PubMed]  

22. S. G. Podorov, K. M. Pavlov, and D. M. Paganin, “A non-iterative reconstruction method for direct and unambiguous coherent diffractive imaging,” Opt. Express 15, 9954–9962 (2007). [CrossRef]   [PubMed]  

23. S. G. Podorov, A. I. Bishop, D. M. Paganin, and K. M. Pavlov, “Mask-assisted deterministic phase–amplitude retrieval from a single far-field intensity diffraction pattern: Two experimental proofs of principle using visible light,” Ultramicroscopy 111, 782–787 (2011). [CrossRef]   [PubMed]  

24. O. Raz, O. Schwartz, D. Austin, A. S. Wyatt, A. Schiavi, O. Smirnova, B. Nadler, I. A. Walmsley, D. Oron, and N. Dudovich, “Vectorial phase retrieval for linear characterization of attosecond pulses,” Phys. Rev. Lett. 107, 133902 (2011). [CrossRef]   [PubMed]  

25. O. Raz, N. Dudovich, and B. Nadler, “Vectorial 1-d signals,” IEEE Trans. Signal. Proces. 61, 1632–1643 (2013). [CrossRef]  

26. C. Song, H. Jiang, A. Mancuso, B. Amirbekian, L. Peng, R. Sun, S. S. Shah, Z. H. Zhou, T. Ishikawa, and J. Miao, “Quantitative imaging of single, unstained viruses with coherent x rays,” Phys. Rev. Lett 101, 158101 (2008). [CrossRef]   [PubMed]  

27. M. A. Pfeifer, G. J. Williams, I. A. Vartanyants, R. Harder, and I. K. Robinson, “Three-dimensional mapping of a deformation field inside a nanocrystal,” Nature 442, 63–66 (2006). [CrossRef]   [PubMed]  

28. J. N. Clark, L. Beitra, G. Xiong, A. Higginbotham, D. M. Fritz, H. T. Lemke, D. Zhu, M. Chollet, G. J. Williams, M. Messerschmidt, B. Abbey, R. J. Harder, A. M. Korsunsky, J. S. Wark, and I. K. Robinson, “Ultrafast three-dimensional imaging of lattice dynamics in individual gold nanocrystals,” Science 341, 56–59 (2013). [CrossRef]   [PubMed]  

29. K. Gaffney and H. Chapman, “Imaging atomic structure and dynamics with ultrafast x-ray scattering,” Science 316, 1444–1448 (2007). [CrossRef]   [PubMed]  

30. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35237 (1972).

31. V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A 20, 40–55 (2003). [CrossRef]  

32. I. McNulty, J. Kirz, C. Jacobsen, E. H. Anderson, M. R. Howells, and D. P. Kern, “High-resolution imaging by fourier transform x-ray holography,” Science 256, 1009–1012 (1992). [CrossRef]   [PubMed]  

33. H. M. L. Faulkner and J. M. Rodenburg, “Movable aperture lensless transmission microscopy: A novel phase retrieval algorithm,” Phys. Rev. Lett. 93, 023903 (2004). [CrossRef]   [PubMed]  

34. J. M. Rodenburg, A. C. Hurst, A. G. Cullis, B. R. Dobson, F. Pfeiffer, O. Bunk, C. David, K. Jefimovs, and I. Johnson, “Hard-x-ray lensless imaging of extended objects,” Phys. Rev. Lett. 98, 034801 (2007). [CrossRef]   [PubMed]  

35. E. J. Candes, Y. C. Eldar, T. Strohmer, and V. Voroninski, “Phase retrieval via matrix completion,” SIAM. J. Imaging. Sci. 6, 199–225 (2013). [CrossRef]  

36. I. Waldspurger, A. DAspermont, and S. Mallat, “Phase recovery, maxcut and complex semidefinit programming,” arXiv.org/1206.0102 (2012).

37. H. Deighton, M. Scivier, and M. A. Fiddy, “Solution of the two-dimensional phase-retrieval problem,” Opt. Lett. 10, 250–251 (1985). [CrossRef]   [PubMed]  

38. J. D. Jackson, Classical Electrodynamics (Wiley, 1975).

39. Y. M. Bruck and L. Sodin, “On the ambiguity of the image reconstruction problem,” Opt. Commun 30, 304–308 (1979). [CrossRef]  

40. F. Hüe, J. Rodenburg, A. Maiden, and P. Midgley, “Extended ptychography in the transmission electron microscope: Possibilities and limitations,” Ultramicroscopy 111, 1117–1123 (2011). [CrossRef]   [PubMed]