Coherent diffractive imaging (CDI) is an experimental technique that measures diffraction intensities and utilizes a numerical phase retrieval algorithm to recover the sample image. These algorithms typically impose constraints on the field estimate using a series of projections applied in an iterative fashion [1]. The method typically relies on a coherent incident field such that far-field propagation is proportional to the Fourier transform. However, short-wavelength sources such as high-order harmonic generation or undulator-based sources have spectra that intrinsically comprise several wavelengths [2,3] or are developed explicitly to have these properties for pump–probe and multi-wavelength experiments [4–6].

Several coherent lensless techniques such as a two-pulse CDI scheme [7], Fourier transform holography [8], coherent modulation imaging [9], and ptychography [10,11] have been developed for polychromatic data. Polychromatic CDI experiments have been demonstrated with the aim of increasing beamline efficiency and reducing exposure times [12–16]. These experiments have demonstrated that reduced temporal coherence can be accounted for by the numerical algorithm. Drawbacks of these methods are the additional algorithmic complexity and computational resources needed for the propagation of polychromatic fields. More importantly, these methods recover a single exit surface wave (ESW) or transmission function, making them inapplicable to general samples.

The proposed method lifts these restrictions on the sample by recovering each ESW separately while retaining the simple experimental setup associated with monochromatic CDI experiments. This is accomplished by exploiting the wavelength-dependent scaling of the real-space coordinate to obtain different supports for each field. The reconstructions are not completely independent, but are coupled through their common support shape, known spectrum, and diffraction measurements. While we describe this method for multi-wavelength CDI, a similar approach could be applied to other situations in which the data comprise an incoherent superposition of intensities, e.g., when the scattered field contains different polarization states [17,18] or time delays beyond the source’s longitudinal coherence. Additionally, the analysis is performed for the 2D case, but is equally applicable to higher dimensions when the real and reciprocal spaces are related through an $ n $-dimensional Fourier transform.

In this Letter, we consider sources whose spectra comprise several distinct wavelengths, ${\lambda _i}$, for $i \in \{1,2 \ldots P\}$, which are sorted such that ${\lambda _1} \lt {\lambda _2} \lt \ldots {\lambda _P}$. The ESWs associated with ${\lambda _i}$ are denoted by ${u_i}$, and we use $i = 0$ as a reference index that can be chosen freely from any of the $P$ indices. The simplest case is treated where each ${u_i}$ is fully coherent with a spatially localized sample that can be treated within the projection approximation. Using the paraxial approximation, the far field is then related to the Fourier transform of the ESWs, ${\tilde u_i}$, by

In this method, the Fourier transform operator, ${\cal F}$, is used for far-field propagation. Taking the inverse Fourier transform of Eq. (1), we find

As we can see from this equation, the spatial coordinate is scaled by the wavelength-dependent factor, ${\varphi _i}$, resulting in reconstructions with different pixel sizes. This results in a wavelength-dependent resizing of the support domains such that

Next, we describe several projections that form the basis for the method. The projections are determined by the known source spectrum, ${S_i}$, and the measured polychromatic far-field diffraction intensities, ${I_m}(y) = I(y) + n(y)$, where $n(y)$ is an unwanted noise term. The supports can be determined throughout the iterative process using a procedure similar to the one described in [19]. Each ${\check{u} _i}$ must fulfill separate support constraints described by Eq. (4). These constraints are fulfilled through the application of the support projection separately to each ${\check{u} _i}$. The support projections at the $k$th iteration are given by

Next, each $\check{u} _i^k$ is normalized to the measured polychromatic amplitude data, $g(y) = \sqrt {{I_m}(y)}$, according to

In Fourier space, the magnitude projection is given by

The projections just described form the basis of the multi-wavelength phase retrieval method. The block diagram in Fig. 1 illustrates the overall structure of the algorithm. It is apparent that the number of unknowns has increased substantially compared to a monochromatic CDI experiment, so next, we discuss under what conditions we can expect to find a unique solution.

 

Fig. 1. Block diagram of multi-wavelength phase retrieval using an alternating projection approach. The algorithm is stopped after a fixed number of iterations and produces $P$ separate reconstructions, ${\check{u} _i}$.

Download Full Size | PDF

Uniqueness is an important aspect of any inverse problem and is discussed through the adaptation of the concepts in [20] to the multi-wavelength situation. By uniqueness, we mean a unique equivalence class consisting of solutions related through translation, inversion, and global phase rotation, $[u] = \{u:u(x) \sim u(x + a),u(x) \sim {u^*}(- x),u(x) \sim {e^{{i\alpha}}}u(x)\}$, where $\alpha \in {\mathbb R}$, $a \in {{\mathbb R}^2}$, and ${u^*}$ denotes the complex conjugate of $u$. For theoretical considerations of uniqueness, see [21–23].

A measure of sufficiency of the constraints is given by the constraint ratio $\Sigma$ [20]. This is a measure of the number of constraints relative to the number of unknowns, which, for the monochromatic case, is given by ${\Sigma _{{\rm mono}}} = \frac{{|A|}}{{2|\Omega |}}$, where $|A|$ and $|\Omega |$ are the size of the autocorrelation and sample supports, respectively. As $\Sigma$ increases, it becomes more likely that a unique solution can be found. In the polychromatic case, $\Sigma$ must include the total size of the autocorrelation function as well as the total number of coefficients within all of the supports. The multi-wavelength constraint ratio becomes

For convex $\Omega$, $|\bigcup\limits_{i = 1}^P {A_i}| = |{A_1}|$, and it is evident that increasing the number of wavelengths results in the reduction of $\Sigma$. This can make the solution more difficult to find or potentially non-unique when $\Sigma \lt 1$. This can be compensated for by using a non-convex or multiply connected support geometry. For a multiply connected support, each ${A_i}$ becomes partially separated with respect to each other, which leads to an increase in $|\bigcup\limits_{i = 1}^P {A_i}|$ and larger $\Sigma$ values. Ultimately, having a large discrepancy between the wavelengths and separated support components will strengthen the support projection in Eq. (5) and increase $\Sigma$, resulting in improved algorithm performance. These concepts are illustrated numerically in the next section through a series of simulations with different support shapes and varying spectra.

Results from several numerical simulations are provided to illustrate the influence that the support, source spectrum, and degree of complexity within the ESWs have on the algorithm performance. The degree of complexity, ${\phi _{{\max }}}$, is the amount of phase within each ESW. As we will see, larger values of ${\phi _{{\max }}}$ typically result in lower reconstruction qualities. Both the support shape and source spectrum affect $\Sigma$, which, as we will see, has a significant impact on the reconstruction qualities.

The images shown in Figs. 2(a) and 2(b) were used to create the two ESWs, ${u_1}$ and ${u_2}$. These $256 \times 256$ images and their 90° rotated versions were used to make the amplitude and phase distributions. The ${\phi _{{\max }}}$ within each ESW is varied from zero to $2\pi$. An example of ESWs with ${\phi _{{\max }}} = 2\pi$ are shown in Figs. 2(c) and 2(d). The complex-valued images are visualized using the colormap shown in Fig. 2(c), which maps the magnitude and phase to value and hue within the image. The four different support shapes are shown by the reconstructions in Fig. 3. The support shapes consisting of one rectangle, two rectangles, one triangle, and two triangles have constraint ratios: ${\Sigma _{{\rm mono}}} = \{2.0,2.8,2.7,4.2\}$, respectively. Experimentally, the support domain is often defined by an opaque sample mask. The ESWs are then far-field propagated to the detector plane; the individual contributions from each ${\hat u_i}$ is shown in Fig. 2(e). The total intensity was calculated using the incoherent sum of intensities according to Eq. (2). Noise was introduced into the data by including Gaussian white noise and quantization of the intensity values. The signal-to-noise ratio was set to 60 dB, calculated according to ${\rm SNR} = 10\log \parallel\! {I_{{\rm exact}}}\!\parallel /\!\parallel\! {I_m} - {I_{{\rm exact}}}\!\parallel$, where ${I_{{\rm exact}}}$ and ${I_m}$ correspond to the ideal and corrupted intensities, respectively. An example of the “measured” polychromatic data is shown in Fig. 2(f).

 

Fig. 2. (a), (b) Paraglider and sailboat images used to simulate the ESWs ${u_1}$ and ${u_2}$, respectively. Red lines indicate the support boundaries. (c), (d) Complex-valued ESWs created using the images from (a), (b). (e) Polychromatic diffraction pattern showing the individual contributions from ${\lambda _1}$ and ${\lambda _2}$ in red and green, respectively, with ${S_i} = \{0.6,0.4\}$. (f) Polychromatic data used by the algorithms.

Download Full Size | PDF

 

Fig. 3. Reconstructions with varying support shapes and degrees of complexity for both wavelengths. The last row shows the exact ESWs with the same supports as (e)–(h). The leftmost (rightmost) two columns correspond to ESWs with ${\phi _{{\max }}} = \pi /2$ (${\phi _{{\max }}} = 2\pi$). The magnitude and phase of the ESWs are mapped to the value and hue within the image using the colormap from Fig. 2(c).

Download Full Size | PDF

The reconstruction procedure is outlined in Fig. 1. We used an algorithm that alternated between 50 iterations of the error reduction method and 100 iterations of a hybrid input–output method [24] for 1500 iterations followed by 500 iterations of the error reduction method for a total of 2000 iterations. The monochromatic projections were replaced with their multi-wavelength counterparts described in Eqs. (5)–(7). A demo version of the code is available at https://bitbucket.org/malmy002/multiwavelength_cdi. Other methods such as the relaxed average alternating reflections [25] or difference map [26] could have been used as well. Fixed supports were used to isolate the algorithm performance from any support determination scheme. A multi-wavelength version of the shrinkwrap method can be used to determine the supports iteratively using the following procedure: first, determine ${\Omega _0}$ by thresholding a Gaussian smoothed version of $|{u_0}|$, then resize ${\Omega _0}$ according to Eq. (4) to obtain the remaining ${\Omega _i}$.

Reconstructions corresponding to each of the four supports and two degrees of complexity (${\phi _{{\max }}} \in \{\pi /2,,2\pi \}$) are shown in Fig. 3. The spectrum is such that ${\lambda _1}/{\lambda _2} = 0.75$ for all reconstructions. The effect of the composition of ${\varphi _i}(x)$ with ${\check{u} _i}$ in Eq. (3) is clearly visible within the reconstructions, as the ${\check{u} _1}$ reconstructions appear larger compared to ${\check{u} _2}$. The reconstructions with triangular supports [Figs. 3(i) and 3(l)] show a clear improvement over the ones with rectangular supports [Figs. 3(a) and 3(d)]. The difference can be explained using ${\Sigma _{{\rm mono}}}$, which corresponds to 2.0 and 2.7 for the rectangular and triangular supports, respectively. In addition, the error increases for samples with higher degrees of complexity, which can be seen by comparing the reconstructions of the first and last two columns in Fig. 3. These observations are supported by the error plots in Fig. 4 where each point corresponds to the average error of 20 independent reconstructions. The relative error was calculated with respect to the known images using the equation $E = \sum\nolimits_{i = 1}^P \parallel \alpha |{u_i}| - |u_i^{{\rm exact}}|\parallel /\parallel u_i^{{\rm exact}}\parallel$, where $\alpha = \parallel\! u_i^{{\rm exact}}\!\parallel /\parallel\! {u_i}\!\parallel$ is a normalization factor. The absolute values are compared to avoid issues associated with global phase offsets between the two images. Because the supports were several pixels larger than the exact supports, each ${u_i}$ was registered to ${u_{{\rm exact}}}$ before computing the error. This registration procedure is described in [27].

 

Fig. 4. (a) Relative error as a function of ${\lambda _1}/{\lambda _2}$ for different support geometries. The solid and dashed lines correspond to ${\phi _{{\max }}} = 0$ and ${\phi _{{\max }}} = 2\pi$, respectively. (b) Relative errors as a function of the constraint ratio. The solid and dashed lines correspond to ${\phi _{{\max }}} = \pi /2$ and ${\phi _{{\max }}} = 2\pi$, respectively.

Download Full Size | PDF

The plots in Fig. 4(a) show that the error decreases as ${\lambda _2}$ deviates further from ${\lambda _1}$, which can be attributed to the relative increase in $\Sigma$. As we would expect, reconstructions associated with ${\phi _{{\max }}} = 0$ (solid lines) have lower errors than those with ${\phi _{{\max }}} = 2\pi$ (dashed lines). Inspection of the errors associated with ${\phi _{{\max }}} = 2\pi$ for the two triangular supports shows that a minimum is reached at ${\lambda _1}/{\lambda _2} = 0.85$. The reason for this is visible in the reconstruction in Fig. 3(o), which shows horizontal-line artifacts that appear at the boundaries of ${\Omega _2}$. The reason for the increased prominence of this artifact is unclear, but it is important to note that it does not appear in the reconstructions with smaller ${\phi _{{\max }}}$ values. Last, the drop in errors in Fig. 4(b) near ${\Sigma _{{\rm mono}}} \sim 2.7$ stems from the relative increase in $\Sigma$ compared to ${\Sigma _{{\rm mono}}}$ for multiply connected supports.

In this Letter, we have developed a general method that extends phase retrieval for multi-wavelength data. Here, we exploit the wavelength dependence of the support domains to separate and recover each spectral component. In general, the introduction of additional unknowns results in lower constraint ratios placing limitation on the number of wavelengths. The method enables single-shot multi-wavelength CDI while retaining the simple experimental geometry of a monochromatic plane-wave CDI experiment, opening up new possibilities for future CDI experiments.