1. Introduction

Determining the phase of the optical field remains one of the most challenging problems in classical optics. These kinds of problems are ubiquitous in various imaging applications: coherent diffraction imaging (CDI), astronomical imaging, wavefront sensing, optical communications, fluorescence tomography reconstruction among others [1–5]. Common approaches to retrieving the phase use either interferometric methods that involve recording an interferogram with a reference beam or non-interferometric methods wherein iterative algorithms are commonly used to derive the phase information from intensity measurements. For the second category of methods, the Gerchberg-Saxton (G-S) algorithm [6,7] and its off-shoots such as the hybrid-input output (HIO) algorithm and the Error-reduction (ER) algorithm [8–10] are well-known approaches.

Recently, machine learning approaches have also been widely adopted to retrieve phase information from intensity-only measurements [11,12]. Despite the success of these approaches, they do require large amounts of data to train a neural network and they can be susceptible to hallucinations i.e. the appearance of features in the final image that can be loosely understood as artefacts due to a mismatch between the training dataset and the actual image. For this reason, there is still interest in the development of robust inverse retrieval methods.

Although the main algorithms offer minimal computational complexity, they aim at solving a nonlinear and a highly ill-posed inverse problem. There are inherent ambiguities associated with phase retrieval that prevent the iterative methods from converging to the correct solution. Non-uniqueness arises from the fact that shifted phases, added global phase and mirrored phases can give multiple solutions to a given phase reconstruction problem [13,14].

In particular, the complex-valued function $g(x,y)$ and its mirrored conjugate function ${g}^\ast (-x,-y)$ both have the same Fourier magnitude, and can be reconstructed with equal probability in any run of HIO algorithm initiated by a random phase guess. In an another possible outcome, the reconstructed image shows features of both the ground truth $g(x,y)$ and its twin-image ${g}^\ast (-x,-y)$. This is referred to as the stagnated mode, it is not trivial to escape it and it most commonly occurs when the object has centrosymmetric support [14]. In one of the methods to overcome this, truncating the object support to a non-centrosymmetric window is employed for the first few iterations. This helps one of the two solutions to dominate. The algorithm subsequently switches to the full support which allows convergence to either the upright or the twin image [13]. Alternatively, the reduction of the total variation for the complex-valued function $g(x,y)$ in a few sub-iterations of the standard HIO algorithm can also help overcome stagnation and yield either of the twin images [15,16]. In most of the existing literature however, reconstruction of either the upright $g(x,y)$ or the twin image ${g}^\ast (-x,-y)$ is accepted as the solution of the phase retrieval problem and the problem of retrieval of the actually correct image is not addressed.

Moving away from these approaches where phase recovery algorithms use single shot Fourier intensity measurements, algorithms based on multiple intensity measurements help curb the issue of non-uniqueness up to a large extent. It has been shown that a series of overlapping diffraction patterns measured by scanning an aperture over the object wave function can lead to overcoming the ambiguity in the phase retrieval. [17]. The diverse measurements could also be for example defocussed intensities or an added intensity measurement taken with a spiral phase placed in the path of the target phase [18,19]. However, the requirement for multiple measurements can introduce experimental complexities or may prove infeasible when the object changes its configuration before the different measurements can be made. Guo et al. proposed perturbing the phase guess in the object domain at each iteration, with a random spatial phase map whose value decreases linearly to escape the stagnation [20].

It has been known that the non-uniqueness is more pronounced for one-dimensional objects [13]. Numerical and theoretical studies show that the Gerchberg-Saxton method applied to non-classical state of light measurements (i.e. photon correlation measurements) can yield unique reconstructions for one-dimensional object which is attributed to a larger Hilbert space that well constrains the algorithm [21]. However, photon correlation measurements for higher number of photons/modes can make the practical realization of this method challenging.

In this work, we address the problem of removing the twin-image ambiguity for phase recovery problems based on single-shot Fourier intensity measurements. We demonstrate that the mathematical properties of the stagnated solution obtained in any run of a typical HIO algorithm can be leveraged to remove the stagnation.

2. Superposed phase from stagnated image

We assume the complex transmittance of an optically thin object is given by:

Here, $A(x,y)$ is the amplitude (or modulus) response and $\phi (x,y)$ is the phase response. If this object is illuminated by a coherent field denoted by $U(x,y)$, then the measurement $O(f_x,f_y)$ registered by a detector placed in the Fourier plane is given by:

Here $\mathscr {F}[.]$ denotes the Fourier transform operation. For our analysis, we assume weakly absorbing objects, hence $A(x,y)$ is unity. The illumination is considered to be normally incident uniform illumination, such that $U(x,y) = 1$. Under these approximations then, and denoting the Fourier magnitude as $|G(f_x,f_y)|$, we can write:

The phase retrieval problem involves finding the phase function $\phi (x,y)$, given the measurement $O(f_x,f_y)$. The G-S/HIO and related algorithms aim to solve this iteratively by propagating the guess phase back and forth between the Fourier and the image domain where we alternately enforce the magnitude constraint and any a priori constraints in the two respective domains.

In any run of G-S/HIO algorithm, it often happens that the solution may stagnate at a linear combination of the function $g(x,y)$ and its twin ${g}^\ast (-x,-y)$. The algorithm fails to move the solution beyond the stagnation point for further iterations and hinders the reconstruction of the correct phase.

However, these stagnated images show an interesting feature: the phase corresponding to the stagnated function when rotated by $180^{\circ }$ (i.e. subjected to point reflection) and subtracted from itself gives a function which is equal to the linear combination of the ground truth phase and its point-reflected version. We can show this analytically. Based on the phenomenological observation that the stagnated image is composed of a superposition of the image and its twin (and that both are identifiable), we can then write the stagnated complex valued function $s(x,y)$ in the near field as a convex combination of the two:

The phase $\alpha (x,y)$ corresponding to this stagnated field is given by:

The phase corresponding to the rotated stagnated field can then be represented as:

Subtracting these two phases, followed by algebraic manipulations, gives:

We point out that the above analysis is independent of the value of the weight $t$ used to represent the stagnated image in Eq. (4), implying that the degree to which the twin image is present does not affect the evaluation of the superposed phase. Our proposal is that this superposed phase can be used as an added information to guide the convergence of the iterative phase retrieval algorithm towards a meaningful solution.

As a side note, we point out that Eq. (4) is in keeping with Ref. [14] where the authors pointed out that the Fourier domain divides into two types of regions: in some, the phase $\alpha (x,y)$ associated with $G(f_x ,f_y )$ (and hence $g(x,y)$ in the object domain) is approximately retrieved and in the second type, the phase $-\alpha (x,y)$ associated with $G^* (f_x ,f_y )$ (and hence corresponding to $g^* (-x,-y)$ in the object domain) is retrieved. Here we then express this recovered solution in the Fourier domain as:

Here $M_1 (f_x,f_y )$ and $M_2 (f_x,f_y )$ are the two binary masks that have the value 1 in the domain where the upright image and the flipped image is present, respectively (and hence complement of each other). Inverse Fourier transforming this to get the expression for the stagnated field in the object domain gives:

Here $*$ represents convolution. Reference [14] shows that in all the examples considered, these two regions, $M_1 (f_x,f_y )$ and $M_2 (f_x,f_y )$ largely form two well-separated, and large areas in the Fourier domain. Then, their inverse Fourier transforms $m_1 (x,y)$ and $m_2 (x,y)$ will be small, i.e. to first approximation, $\delta$ functions that convolve $g(x,y)$ and $g^*(-x,-y)$, respectively. The stagnated phase image will therefore appear to be a linear superposition of the flipped and the upright images. This is indeed, exactly what one phenomenologically observes. This therefore leads us to our phenomenological relation for the stagnated image given by Eq. (4). The statistical validity of Eq. (4) is studied over different random trials and shown in Fig. 1(f).

 

Fig. 1. Illustration of phase stagnation. (a) Test grayscale phase image $\phi (x,y)$. (b) Fourier magnitude for the test phase image, plotted in log scale. (c) Stagnated phase obtained after 800 iterations of the standard HIO algorithm. (d) Ground truth for the superposed phase $\phi (x,y)-\phi (-x,-y)$. (e) Superposed phase evaluated by subtracting the point inverted phase of the stagnated field from itself. (f) Relative error between the superposed phase obtained from the stagnated image and the ground truth for 100 different random initiations of the HIO algorithm. The maximum and the minimum relative error values over 100 trials are 0.21 and 0.07 respectively, with a mean of 0.14, showing good fidelity with the ground truth and statistical stability over many repetitions with different random initializations

Download Full Size | PDF

3. Augmented iterative phase retrieval (AIPR)

The proposed method comprises of a certain number (e.g. 1400) iterations out of which a first set (e.g. 800 iterations) follow the standard HIO update where a random initial guess, the known Fourier magnitude $|G(f_x,f_y)|$ (alternatively, $\sqrt {O(f_x,f_y)}$) and the object support are fed as inputs. We purposely choose our test phases to have centrosymmetric support as these kind of images are known to be more prone to stagnation [15]. During this stage the algorithm explores the solution space and yields either the stagnated field or the correct solution. At the end of 800 iterations, the phase of the output is subjected to point inversion and subtracted from itself to get the superposed phase. The phase corresponding to the 800th update is reserved as the input for the next iteration steps, where augmentation starts.

Figure 1 demonstrates the generation of the superposed phase from the stagnated function for a grayscale test phase image $\phi (x,y)$. Fig. 1(c) shows the stagnated phase reconstructed from a run of 800 iterations of HIO algorithm. We test the validity of the equality in Eq. (7) by running 100 trials of such 800 HIO iterations with different random initiations and evaluating the superposed phase from the output. The relative error between thus generated superposed phase and the ground truth superposed phase (shown in Fig. 1(d)) for each of the 100 trials is shown in Fig. 1(f). The consistently low values indicate that the superimposed phase can be produced with high fidelity from the stagnant field functions, matching well with its respective ground truth. We point out that the 100 trials here are only shown to verify the statistical robustness of the superposed phase generation method.

Our algorithm is based on the finding that the stagnated mode, post the 800 iterations can be overcome by perturbing the phase guess obtained henceforth after each HIO update, with a combination of two superposed phases. In line with the description above, these are obtained by applying the operators $\mathscr {R}_-$ that takes as its input, any 2-D phase function $\beta (x,y)$, rotates it by $180^{\circ }$ and subtracts the rotated input phase function from the original input:

We call the update steps as the augmentation steps. As is the case with HIO algorithm, the feedback parameter lies between $[0,1]$ [9] and it was empirically fixed to 0.8 for our analysis. The detailed steps for the $(n+1)th$ iteration of the algorithm can be described via the pseudocode described in Algorithm 1.

Algorithm 1. Augmented iterative phase retrieval

View Table

4. Simulation results

4.1 Effect of the augmentation steps

In order to understand the impact of the augmentation steps, we illustrate the alterations they introduce in the trajectory of the algorithm. To make sure that the Nyquist sampling is followed for all the intensity images considered in our analysis, we have zero padded our images such that the computational window is of the size of 512 $\times$ 512 pixels and the object support size covers the central 250 $\times$ 250 pixels. The results presented henceforth show the central 250 $\times$ 250 pixels region. The algorithm presented here aims to retrieve complex valued objects and the test images shown henceforth are used as phase of the input field. We use the Fourier magnitude error and the phase error for our analysis which are defined, respectively, as:

Here $M$ is the total number of pixels and $\hat {\phi }_{n}$ denotes the phase guess at the $n$th iteration number. The blue curve in Fig. 2(a) shows the error metric $E_{1}$ for the standard HIO algorithm which has stagnated in a local minimum. The inset shows the stagnated phase image at the 800th iteration. Application of the augmentation steps from this point onward brings a substantial change in the error value. The red curve shows the trajectory for the case when augmentation is stopped after the 820th iteration (hence 20 augmentation steps) and subsequently HIO is continued. We observe that the error starts to proceed towards a lower value (red curve). However, if we continue to apply the augmentation steps, it leads to a larger error value. In this range further implementation of augmentation steps do not bring substantial change in the phase guess, eventually leading to a wrong solution. Figure 2(b) shows the detailed plot of the changes introduced due to seventy augmentation steps post the 800 HIO iterations. We can broadly divide the course of evolution in three regions. We see that there is a small change in the error metric for the first few iterations (represented by Region I in the figure). The phase error plots in Fig. 2(c-d) show that the repeated small changes induced by the augmentation can eventually bring a substantial change in the updates, pushing the phase guess out of the local minimum and into Region II where we see a sudden large increase in the error. This sudden increase after several AIPR iterations is a universal feature that we observed in all tested image retrievals. Continuing with the AIPR augmentation will then continue to increase the error until this settles to some roughly constant value (Region III) and consequently provides wrong images. However, once the solution has been pushed out of the local minimum, returning in Region II to the standard HIO algorithm brings convergence towards a new solution with a reduced error and which we found in all cases to lead to the correct solution.

 

Fig. 2. Comparative error evolution with same random initialization. (a) Fourier magnitude error $E_{1}$ for three cases: standard HIO algorithm, AIPR algorithm with augmentation applied after step 800 but only for limited steps and then followed by a return to HIO, and AIPR algorithm with augmentation applied continuously after step 800. The inset shows the reconstructed phase image at different stages of the three cases. (b) Detailed section from error plot of (a) showing the first seventy augmentation steps (800-870) that highlight the three regions I, II and III. (c) Corresponding phase error for the three cases. (d) Detailed section from error plot of (c) showing the first seventy augmentation steps (800-870).

Download Full Size | PDF

4.2 Reconstruction for test phase images

We test the performance of the proposed algorithm via numerical simulations on gray-scale images taken as test phase. The phase values are scaled in the range of $[0,\pi /2]$. We implement 1400 iterations of the proposed augmented iterative phase retrieval algorithm wherein the superposed phase was constructed after the first 800 iterations. The augmentation steps were implemented at this point for a number of steps before switching back to the standard HIO algorithm again. The number of augmentation steps $N_{aug}$ depends on the phase image and the initial guess and it varied from 17 to 19 for the four cases considered here.

Figure 3 shows the reconstructions of four test phase images from their Fourier magnitudes using AIPR, compared with the reconstructions using the standard HIO algorithm. For comparing each of the reconstructed images with the ground truth phase function, we use the Pearson correlation co-efficient, which for two generic functions $a$ and $b$, is defined as:

 

Fig. 3. Reconstruction results for different test phase images. (a-d) Grayscale test phase images. (e-h) Stagnated reconstructions after 800 iterations of the standard HIO algorithm. The PCC values are respectively, 0.6213, 0.4877, 0.7795, and 0.6491. (i-l) AIPR algorithm reconstructions carried out with same random initiations as that used for standard HIO algorithm reconstructions in the second column. PCC values are respectively, 0.9610, 0.9998, 0.9958, and 0.9873. (m-p) Fourier magnitude error $E_{1}$ for the four cases.

Download Full Size | PDF

We point out that this approach can lead to the unambiguous upright solution if we construct a combination of superposed phases obtained from a) the final output phase and b) from the stagnated phase. To this end, in any run of the proposed method the final output phase, e.g. $\phi _{1400}$ (which could be either upright or flipped) can be subjected to the operator $\mathscr {R_+}$. For both the cases of the flipped and the upright image, the result will be the same, i.e. $\phi (x,y)+\phi (-x,-y)$. As mentioned before, the operator $\mathscr {R_-}$ applied on the phase output of the 800th iteration gives $\phi (x,y)-\phi (-x,-y)$. It can be seen that a final combination of these two will yield the upright solution. Specifically, the upright image can be obtained as:

The method as it stands has not been optimized as to get the precise number of the augmentation steps. However, as demonstrated in Fig. 2, it was established phenomenologically that the first few steps did not bring the required change and moving beyond 20-25 augmentation steps yielded another local minimum, with larger Fourier magnitude error. One needs to be in the Region II where we see considerable change in the phase guess as compared to the previous guess. Once in region II, the decision probably could be automated by running the HIO for a few steps and then evaluate the degree of symmetry of the phase image. This forms the subject of future studies based on this method.

4.3 Reconstruction for larger phase range images

The range chosen for the analysis so far was motivated by the fact that the $tan^{-1}$ function in Eq. (7) holds when the test phase image lies in the range $[0,\pi /2)$, or equivalently the superposed phase lies in the range $(-\pi /2,\pi /2)$. There will be a larger uncertainty in superposed phase evaluation for wider phase range, owing to phase wrapping. To minimize this error, we do a modulo operation on all those values of the evaluated superposed phase that are outside the range $(-\pi /2-\delta \epsilon,\pi /2+\delta \epsilon )$. The $\delta \epsilon$ here is introduced as tolerance, keeping in mind that this being a numerical process there will be some error associated with each pixel value of the superposed phase and a hard boundary of $(-\pi /2,\pi /2)$ for performing modulo operation can alter the otherwise approximately correct values at the limits of this boundary.

Figure 4 shows the results for the test phase image in $[0,4\pi /5]$ range. Figure 4(e) shows the relative error between the superposed phase obtained from the stagnated image and the ground truth for 100 different random initiations of the HIO algorithm. As expected, the mean relative error (0.38) is larger than that reported for the case when the test phase image was in the range $[0,\pi /2]$. AIPR was able to reconstruct the right phase image as shown in Fig. 4(f) with a PCC of 0.9705 and the number of augmentation steps being 11.

 

Fig. 4. Results with larger phase range. (a) Test phase image in the range $[0,4\pi /5]$.(b) Stagnated phase obtained after 800 iterations of the standard HIO algorithm.(c) Ground truth superposed phase. (d) Superposed phase evaluated by subtracting the point inverted phase of the stagnated field from itself. (e) Relative error between the superposed phase obtained from the stagnated image and the ground truth for 100 different random initiations of the HIO algorithm. (f) AIPR reconstruction.

Download Full Size | PDF

4.4 Reconstruction for noisy data

We study the effect of noise on the AIPR algorithm by generating Fourier intensity data limited by Poissonian noise for two different values of average photons per pixel. Figure 5(b,f) show the stagnated phases obtained from the noisy Fourier magnitudes for the cases of $10^4$ and $10^3$ photons per pixel, respectively. The phase range used for this study was $[0,4\pi /5]$. The percentage of pixels with zero photoncounts was $24{\%}$ for the higher light level and $48{\%}$ for the lower light level condition. Figure 5(c,g) show the superposed phase generated during AIPR with a relative error of 0.36 and 0.39 respectively for the two cases. AIPR could reconstruct the right phase image with PCC value of 0.94 and 0.91 for the case of $10^4$ and $10^3$ photons per pixel, respectively, although the reconstruction is noisy and there is a faint trace of the twin image in the lower light condition. The number of augmentation steps for the two cases was 14 and 12.

 

Fig. 5. Results with noisy Fourier intensity (a,e) Noisy Fourier magnitudes $G(f_x,f_y)$ (log scale), (b,f) Stagnated phases obtained after 800 iterations of the standard HIO algorithm,(c,g) Superposed phases evaluated during AIPR, (d,h) AIPR reconstructions for average light levels of $10^4$ and $10^3$ photons/pixel, respectively.

Download Full Size | PDF

5. Conclusion

Phase recovery of a complex-valued object from single-shot Fourier magnitude data is an ill-posed problem. We demonstrate that the stagnated near field solutions have mathematical properties which can be used to construct the superposition of ground truth phases. This provides a prior information about the phase object that can be used as a constraint. We show through numerical simulations that the superposed phase can be used in the iterative phase retrieval algorithms to guide the convergence towards the right solution. Notably, this approach not only resolves the stagnation, but also delivers the correct image, i.e. the actual image without any ambiguity with the inverted image solution as is often the case with other approaches. Analysis with noisy data simulate real world data situations and show that the algorithm is still able to come out of the stagnation. The proposed method has its current limitation as lacking an algorithmic way to find the precise number of augmentation steps beyond checking that we are in the Region II of the error evolution. Future implementations of this approach can study algorithmic formulation which have alternations between HIO and augmentations, in line with the algorithms which switch between HIO and the ER algorithms.