Two noise-robust axial scanning multi-image phase retrieval algorithms based on Pauta criterion and smoothness constraint

Cheng Shen; Xuejing Bao; Jiubin Tan; Shutian Liu; Zhengjun Liu

doi:10.1364/OE.25.016235

1. Introduction

Iterative phase retrieval methods play an important role in coherent diffractive imaging (CDI) technique, a useful lensless imaging modality with the potential to achieve resolutions that are not limited by the lens devices. It has undergone a rapid development over the last decade with application in the experiments involving X-rays [1, 2], electrons [3] optical laser [4] and high-harmonic light source [5], and fields like materials science [6], metrology [7] and biology [8–10]. Although these significant advances have been made, the precision of CDI technique has been impaired by the experimental noise since the technique was born [11]. It is a significant task in CDI to ameliorate the noise effect by improving the phase retrieval algorithms.

At the origin of iterative phase retrieval algorithms, Gerchberg-Saxton (GS) algorithm reconstructs the object phase function with a pair of known amplitude distribution at the object plane and at the measuring plane [12]. But it is limited by the low convergence speed and the high sensitivity to initialization as well as the demanding priori knowledge of the object amplitude distribution. Thus, the hybrid input output (HIO) algorithm [13] and the continuous version of HIO (CHIO) [14] were put forward by imposing the support constraint and introducing the feedback control into the iteration. However, the requirement of a tight support and the selection of proper feedback coefficient are still inconvenient. Multiple image phase retrieval algorithms were proposed [15–20] afterwards. They introduced variable optical system parameters to generate multiple measurements and can be categorized into two groups: lateral and axial scanning strategies. Ptychography [1, 2, 18–24] is the representative of lateral scanning techniques and has made lots of achievements theoretically and experimentally. Its basic idea is to create a series of overlapped illumination regions at the sample plane and acquire a redundant data set so that phase retrieval algorithms can reconstruct the phase information. Now ptychography has developed into a relatively independent and mature methodology. Its behavior under noise conditions has been fully discussed in [25–28].

For the axial scanning strategies, the single-beam multiple-intensity reconstruction technique (SBMIR) algorithm [29] and the multi-stage algorithm [30] both sequentially record diffraction or scattering data and process them serially in the iteration update. By contrast, the amplitude-phase retrieval (APR) algorithm [31] deals with the measured intensity patterns in a parallel way by virtue of averaging estimations in the real space. Notable is that these axial scanning multi-image phase retrieval algorithms are not specially designed to reduce the noise effect and the relevant noise study is insufficient. Till now, iterative phase retrieval algorithms designed for the noise robustness include a combination of hybrid input-output (HIO) and error-reduction (ER) algorithm [13], difference map [32], relaxed averaged alternating reflectors [33], noise robust (NR)-HIO [34], and oversampling smoothness (OSS) [35, 36], which are all single-shot phase retrieval algorithms. We will discuss the noise performance of multi-image methods based on APR algorithm under different noise models and in experiment. It needs to be declared that the multi-image phase retrieval methods discussed in this paper exclude ptychography.

In the previous work, the scheme combining single-lens imaging and APR algorithm was verified to hold several advantages, such as simple experimental implementation and easy operation [37]. Based on this system, the APR algorithm and two proposed algorithms are investigated under different kinds of ‘noise’, which contain not only the common experimental noise models but also problems like quantization error and intensity center missing. The missing of intensity data at the center of diffraction patterns is to some extent inevitable due to the beamstop in X-ray experiments and the overexposure in optical experiments. A beamstop is commonly used in X-ray CDI experiments to stop the intense primary beam that has not been diffracted by the sample and collect more high-angle diffraction signals [38]. But it causes the low-frequency signals missed, which can make CDI reconstruction unstable or even fail [39]. It parallels the situation in optical experiments that the limited bit depth of detectors is not able to efficiently capture the large dynamic range of the intensity distribution [40, 41]. The overexposure part corresponding to the low-frequency signals seems like being blocked by a beamstop. Fortunately, what size of missing data will lead to the permanent loss of some structural information on the specimen has been investigated [42]. In this paper, the problem is taken as a general kind of ‘noise’ and analyzed in the algorithm performance.

This paper is organized as follow. First, the improvement scheme is explained based on the mathematical analysis. Then, three kinds of common experimental noises are modeled for the following simulation. Next, three existing algorithms and two proposed algorithms are tested under different simulated noise conditions, including the missing central intensity. It is verified by simulations and a demonstrative experiment that the proposed algorithms hold a consistent superiority over the old ones both in reconstruction fidelity and convergence speed. The work is inspiring in the light of fact that the resolution of CDI technique is limited by the low signal-to-noise level of diffraction data at high scattering angles. Also, reconstructing fine features in objects like biological specimens from noisy experimental data remains still a challenge.

2. Algorithm analysis

The system is shown in Fig. 1(a). It is a simple analyzer and can be mathematically depicted by the extended fractional Fourier transform (eFrFT) [37]. In this paper, we focus on improving the noise robustness of APR algorithm with the Pauta criterion and the smoothness constraint respectively.

Fig. 1 (a) The eFrFT system for phase retrieval. Two isomers of the noise-robust APR algorithm: (b) performing the average operator first, (c) performing the 3σ constraint first.

Download Full Size | PDF

Support is a region to give the estimated size of object and is usually loose since the prior knowledge is hard. For an iterative algorithm, the region outside of support in the final reconstruction should be zero-valued without noise in the measured intensity images. However, during the iteration, it is not always zero-valued, caused by the phase estimate error in the reciprocal space. The phenomenon can be exploited to avoid the stagnation issue. When the noise exists, the non-zero values outside of support during the iteration are a mixture result of the phase estimate error and the amplitude error in the reciprocal space. The latter stems from the measurement noise. It is desirable to exploit the former to avoid stagnation and reduce suffering from the latter simultaneously. Thus, how to distinguish the two different error sources becomes a problem. In a previous work [34], the Pauta criterion is depended on to address the problem. It is shown as

ψ^{(n + 1)} = {\begin{cases} ψ_{\mod}^{(n)} ψ_{\mod}^{(n)} \in S, \\ ψ^{(n)} - β ψ_{\mod}^{(n)} ψ_{\mod}^{(n)} \notin S \cap | ψ_{\mod}^{(n)} | > 3 σ_{noise}, \\ 0 ψ_{\mod}^{(n)} \notin S \cap | ψ_{\mod}^{(n)} | \leq 3 σ_{noise}, \end{cases}

where

ψ^{(n)}

is the iterate at the nth calculation in the real space,

ψ_{\mod}^{(n)}

is a process variable after applying the amplitude constraint to

ψ^{(n)}

and

ψ^{(n + 1)}

is the updated iterate. The symbol

S

represents the support region. The coefficient

β

is taken at 0.9 empirically. If the pixel value outside of support is less than

3 σ_{noise}

, it will be taken as the error brought by the noise and reset to zero. Otherwise, it is deemed to be brought by the phase estimate error during the iteration, which can be exploited to avoid stagnation. According to [34], the standard deviation of the non-zero pixel values outside of support due to noise can be approximated as

σ_{noise} = 0.5

for a mean signal greater than 0.2 photons per pixel in the measured intensity images.

For another hand, the correlation information among the pixels outside of support in the real space can be exploited, which is smoothness constraint [35]. It applies adjustable spatial frequency filters to pixels outside the support at different stages of iteration and searches for a global minimum in the solution space. The protocol is

ψ^{(n + 1)} = {\begin{cases} ψ_{\mod}^{(n)} ψ_{\mod}^{(n)} \in S \cap ψ_{\mod}^{(n)} \geq 0, \\ ψ^{(n)} - β ψ_{\mod}^{(n)} ψ_{\mod}^{(n)} \in S \cap ψ_{\mod}^{(n)} < 0, \\ W (ψ^{(n)} - β ψ_{\mod}^{(n)}) ψ_{\mod}^{(n)} \notin S, \end{cases}

where W represents an adjustable Gaussian filtering operator in the reciprocal space.

When the Pauta criterion is applied to APR algorithm, the NRAPR algorithm, short for noise-robust APR, has two ‘isomers’ due to the order of the average operator in APR algorithm and the 3σ constraint in the Pauta criterion framework, shown in Figs. 1(b)-1(c). The modification of APR algorithm with smoothness constraint, SAPR short for smoothness APR, also comes up against the same problem. Thus the optimum choice needs to be determined. The relevant analysis will be given at the end of Section 3.

3. Noise model

In the CDI experiment, there are three main noise sources. Firstly, shot noise occurs in the photoelectric conversion, associated with the particle nature of light and the discrete nature of electric charge. Secondly, the photoelectric devices, like charge coupled devices (CCD), suffer from a dark current. The pattern of different dark currents can result in a fixed-pattern noise. Dark frame subtraction can remove an estimate of the mean fixed pattern, but there still remains a temporal noise, because the dark current itself has a shot noise. Thirdly, the noise could result from scattering of the substrate, incoherently illuminated sections of the sample, slit scatter or air scatter.

According to their statistics, three corresponding models are established. For the shot noise, typical of multiplicative noise, its simulation procedure is as follows: (i) calculate an ideal measurement assuming a total of photons for the entire pattern; (ii) generate 2-D Poisson distribution with a varying mean equal to each pixel value in the ideal measurement. Then, a sample of the 2-D Poisson distribution can be taken as a pattern corrupted by shot noise. The ideal intensity measurement is symbolized by $I$ and the noisy pattern $P$ . For clarity, we define a normalized difference function

de (A, B) = \frac{\sum_{j = 0}^{N - 1} \sum_{k = 0}^{N - 1} | | A_{j, k} | - | B_{j, k} | |}{\sum_{j = 0}^{N - 1} \sum_{k = 0}^{N - 1} | A_{j, k} |}

to evaluate how the two-dimensional matrix

B_{j, k}

deviates from the truth

A_{j, k}

. So an index

R_{noise} = de(\sqrt{I}, \sqrt{P})

can be introduced to represent the noise level. Obviously, R_noise quantifies the noise level in a contrary way than the signal-to-noise ratio (SNR) does. With the image [43] in Fig. 2(a) as the object, the relationship between shot noise level and the photon number is calculated. As is shown in Fig. 2(b), R_noise plummets when the photon number increases. In the X-ray experiments, the X-ray dose cannot be too heavy considering its devastation to vulnerable biological specimens. Also, the illuminance in optical experiments could be intentionally tuned down to avoid the center overexposure. Thus, it is usual that the photon number is relatively low. A pair of log-scaled ideal pattern

I

and noisy pattern

P

is displayed in Figs. 2(c)-2(d) when the photon number is 10⁹.

Fig. 2 (a) 92 $\times$ 92 pixel tested image, placed in the center of 256 $\times$ 256 pixel background in the calculation. (b) Shot noise level as a function of the photon number. (c) Logarithm of a simulated noise-free diffraction pattern when the photon number is 10⁹. (d) The shot-noise corrupted counterpart of (c).

Download Full Size | PDF

The other two common noises, dark frame subtraction and alien scattering, are both additive noises. To mimic the dark frame subtraction [44], a Gaussian random distribution $D$ is generated and then subtracted from the ideal measurement $I$ . With a variance of 1, the mean of Gaussian random distribution is taken as a varying percentage of $I$ . Accordingly, the index R_noise is calculated as $de(\sqrt{I}, \sqrt{I - D})$ . For alien scattering, it is supposed to be spatially uniform. So a 2-D Poisson distribution $S$ with a mean of between 0.0001 and 1 photon per pixel is generated and then added to $I$ , corresponding to $R_{noise} = de(\sqrt{I}, \sqrt{I + S})$ .

Now, the aforementioned ‘isomer’ issue is analyzed. To evaluate the accuracy and consistency of final reconstructions, the index $E_{r e a l} = de(A, ψ^{(n)})$ is firstly defined, where $A$ denotes the image in Fig. 2(a). Then, a set of noisy diffraction patterns is generated with R_noise approximately 30% when the photon number is 6.79 × 10⁵, which represents a highly noisy condition. Next, they are input into the flowchart (b) and (c) in Fig. 1. E_real after 1000 iterations is calculated for two isomers and shown in Table 1. It can be seen from the value of E_real that the structure (c) outperforms (b) for NRAPR. Thus, the structure (c) is adopted in the following discussion. Similarly, the structure (b) of SAPR is the optimum choice for further discussion.

Table 1. Reconstruction error after 1000 iterations for two isomers of proposed algorithms

View Table | View all tables in this article

4. Results and discussion

In this section, the noise performance under different conditions is investigated for old and proposed algorithms, including NRHIO, OSS, APR, SAPR and NRAPR. Here NRHIO and OSS are both single-shot noise-robust phase retrieval algorithms. They are listed here in order to compare the noise performance between single-shot and multi-image methods.

The test image is Fig. 2(a). In the simulation, the working wavelength is chosen as 632.8nm and the focal length is 200mm. The physical side length of square test images is fixed at 5mm. For multiple measurements, the object distance is taken as 50mm, 70mm, 90mm, 110mm and 130mm sequentially. For single measurement, the object distance is 90mm. Meanwhile, the distance between the object plane and the measuring plane is kept as twice the focal length. It has been testified that the system parameter setting satisfies the sampling theorem.

As for the algorithm parameter setting, the support size of NRAPR and SAPR is 100 × 100 pixels. The coefficient $σ_{noise}$ is kept as 0.5 for NRAPR since the simulated intensity patterns all have a mean signal greater than 0.2 photons per pixel. Also, the best retrieved result from different random initializations is picked up to display in each test. Thus, all simulations follow the Monte Carlo method. Besides, another commonly used error metric $E_{r e c i p r o c a l} = de(M, Ψ^{(n)})$ is calculated, where $Ψ^{(n)}$ is the eFrFT of $ψ^{(n)}$ and $M$ noisy patterns. Therefore, E_reciprocal is an indicator how close the retrieved result is to the noisy measurements in the reciprocal space but not simply the counterpart of E_real.

4.1 Photon number

In this subsection, the shot noise is introduced with R_noise ranging from 5% to 30% according to Table 2. Figures 3 and 4 show the results after 5000 iterations. From Fig. 3(a), it can be clearly seen that two proposed algorithms attain more accurate reconstructions than three existing algorithms do, judging from the value of E_real. The dominance becomes more notable when the noise level goes higher. For instance, when R_noise is 30%, APR fails to get a recognizable reconstruction while the retrieved images of NRHIO and OSS are quite blurry and show less subtle details than NRAPR and SAPR’s, shown in Fig. 3(c). In Fig. 3(a), it can be observed that the result of NRHIO keeps the worst whatever the shot noise level is and SAPR has the best performance overall. Interestingly, the metric E_reciprocal does not show the same trend as E_real. Except NRHIO, the error metrics of other algorithms in the reciprocal space are quite close to each other, which apparently contradicts the results displayed in Fig. 3(c). Since E_reciprocal evaluates how close the retrieved result is to the noisy measurements in the reciprocal space, the smaller E_reciprocal does not necessarily represent a closer proximity to the true object function in the real space. It says E_reciprocal is deceptive and cannot indicate the actual retrieval fidelity.

Table 2. The chosen noise conditions with photon flux density descending

View Table | View all tables in this article

Fig. 3 Reconstruction errors of old and proposed algorithms under the shot noise: (a) $E_{r e a l}$ as a function of the noise level, (b) $E_{r e c i p r o c a l}$ as a function of the noise level, (c) final reconstructions when R_noise = 30%.

Download Full Size | PDF

Fig. 4 The convergence curves of five algorithms when $R_{noise}$ is 5%.

Download Full Size | PDF

In Fig. 4, noticeable is that not only do the two proposed algorithms converge to a more accurate solution, but also have a higher convergence speed, which means that the proposed algorithms have advantages in both accuracy and efficiency. Besides, a phenomenon from Fig. 4 is that the OSS curve oscillates for several hundred iterations after the first convergence, which is attributed to it that the tunable spatial frequency filter in the smoothness constraint gets inappropriate characteristics. Fortunately, its self-adjustment makes the curve converge again soon after.

4.2 Dark frame subtraction

As stated in Section 3, the mean of Gaussian random distribution is taken as 0, 10, 20, 30, 50, 70, 90, 95, 99 and 100% of I, which yields R_noise between 2.05% and 91.49%. Then the bias is subtracted from the perfect measurement when the photon number is 3.75 × 10⁷.

Figure 5 gives the quantitative results of five algorithms’ retrieval after 5000 iterations. Generally, it meets the expectation of an increasing error metrics with R_noise growing. Two single-shot phase retrieval algorithms get the worst reconstruction. Especially when the bias level ranges from 30% to 60%, their E_real are 1-2 orders of magnitude larger than the other three’s. Once the bias level exceeds around 67%, SAPR does not work as well as it did before, and at this time NRAPR gets the most satisfactory retrieval. When the measured patterns are highly corrupted by the bias, APR outperforms all other algorithms. From Fig. 5(b), the trend of E_reciprocal roughly goes with E_real at this time.

Fig. 5 Scatter plots of two error metrics against $R_{noise}$ for the different bias levels.

Download Full Size | PDF

Significantly, APR produces the reconstruction almost as NRAPR does, much better than the two single-shot methods. It suggests the average operator in APR is a powerful constraint under the additive Gaussian noise.

4.3 Alien scattering

Varying the mean of Poisson distribution between 0.0001 and 1 photon per pixel gives the diffraction patterns corrupted by alien scattering with R_noise from 0.69% to 29.21%. The total of photon number is 6.79 × 10⁵.

Figure 6 shows E_real and E_reciprocal after 5000 iterations as function of R_noise. Similarly, OSS and NRHIO arrive at the worst results. The SAPR curve soars after R_noise is beyond 16.22%, despite that it gets the best reconstruction in the lower noise level range. In contrast, the NRAPR curve is growing steadily, which approximately overlaps the APR curve and is surpassed by APR at the high noise level. The observation is close to the one in subsection 4.2. From the perspective of E_reciprocal, APR consistently obtains the best reconstruction from noisy diffraction patterns, which contradicts the truth.

Fig. 6 Scatter plots of two error metrics against $R_{noise}$ for the different alien scattering levels.

Download Full Size | PDF

Considering both dark frame subtraction and alien scattering bring in additive noises, the similarity of conclusions in subsection 4.3 and 4.4 is in expectation. Taking the performance from the whole noise levels into account, NRAPR is the preferable algorithm although APR itself holds an inherent insensitivity to additive noises. Also, it can be concluded that E_reciprocal, the conventional metric used to monitor the iteration, may well malfunction in the actual experiments. Thus, we will put forward a more reliable metric in the following subsection 4.6.

4.4 Quantization error

Apart from the multiplicative and additive noises discussed above, the inevitable quantization error, brought by the analog-to-digital devices, is investigated here. Usually, the more the storage bits are, the more expensive the devices are. And there is a limit to storage bits because of the manufacturing precision. Thus, it is essential to discuss how the quantization error influences phase retrieval. The ideal intensity patterns are represented as 64-bit floating point numbers in the simulation and then converted to the lower-bit data formats from 16-bit to 8-bit numbers.

Figure 7 gives the quantitative comparison of five algorithms. NRHIO holds a high E_real all the time. The reconstruction of NRAPR is the best in the whole range and has a slender advantage over APR. Moreover, two scatter plots are illustrated to reveal the similarity between the trend of OSS and SAPR. They both start at a low error metric, then shoot up and end up with the worst reconstructions. It suggests SAPR is quite sensitive to the quantization error.

Fig. 7 The effect of quantization error: (a) intensity patterns in different data bits when the object distance is 90mm, (b) a mixture of bar charts and scatter plots showing $E_{r e a l}$ against data bits of intensity patterns.

Download Full Size | PDF

4.5 Center missing

In this subsection, the algorithm performance is investigated in the absence of central intensities. In view of the fact that the intensity information within each speckle corresponds to a real-space dimension larger than the sample size, phase retrieval from oversampled diffraction intensities is sensitive only to the number of missing speckles rather than the number of missing pixels [42]. The number of missing speckles can be characterized by

η_{t} = \frac{D_{t} - 1}{2 σ_{t}}, t = x or y,

where

D_{t}

is the number of missing pixels,

t

indicates the lateral directions and

σ_{t}

is the oversampling ratio in the real space [45], defined by

σ_{t} = \frac{density region + no-density region}{density region} .

It has been proved that some sample structural information will be permanently lost to not be recovered by iterative algorithms when

η_{x} = η_{y} = η > 1

[42]. Thus, five algorithms are tested in the condition that

η < 1

and

η > 1

respectively.

Figure 8 provides an overview of the retrieved results. When $η < 1$ , two proposed algorithms have a reconstruction with higher fidelity than others, although the reconstructed images of APR and OSS are acceptable. When $η > 1$ , reconstructions from all iterative algorithms are intensively impaired by missing central intensity, with the appearance of visual artifacts. But the reconstructed images of two proposed algorithms are still the best with less density oscillation. To sum up, two proposed algorithms are still the optimum choice when the missing central intensity exists.

Fig. 8 Reconstructions of old and proposed algorithms when center intensities are missing.

Download Full Size | PDF

4.6 Experimental demonstration

To verify the effectiveness of proposed algorithms, the optical phase retrieval experiments are conducted. The experimental setting is shown in Fig. 9(a). A collimated linearly polarized 532nm laser beam (MW-SGX, Changchun Laser Optoelectronics Technology) is employed to illuminate the aperture mask (sample) and then focused by a lens with the effective focal length 100mm, followed by a 2016 × 2016 CCD (GS3-U3-41S4M, Point Grey Research). The distance between the sample and the lens is 40.81mm. We install the CCD on a precision translation stage (M-403, Physik Instrumente) driven by a servo controller (Mercury C-863, Physik Instrumente). Then, several images are recorded at consecutive distances, which are 70mm, 80mm and 90mm away from the lens. To provide the ground truth for visual comparison, the sample is imaged on an automatic computerized numerical control (CNC) measurement system (Smart Scope MVP 200, Optical Gaging Products) after the experiment. As is shown in Fig. 9(b), it is back illuminated with LED profile light (collimated, green) and the white scale bar corresponds to 0.84mm.

Fig. 9 Optical CDI experiments: (a) the experimental setup for multi-image phase retrieval, (b) a sample image acquired from an automatic CNC measurement system after the experiment, (c) measured intensity patterns at the imaging distance of 70mm, 80mm and 90mm respectively, cropped to 800 × 800 pixels, (d) reconstructed images from the old and proposed algorithms, cropped to 650 × 650 pixels for visibility.

Download Full Size | PDF

In the experiment, the measured intensity patterns are contaminated by the mixed noise. The laser is tuned down to create the dark-field illumination, bringing in the shot noise. To give a rough estimate of the shot noise level, a dual-channel light power meter (2832-C, Newport Corporation) is used to gauge the instantaneous light power P. By setting the same camera gain and recording the exposure time T, the photon number corresponding to a camera image with some specific gray value sum is calculated as $P \times T$ divided by the light quantum. Then the photon number corresponding to a noisy intensity pattern can be estimated proportionally according to its gray value sum. All three measured noisy patterns are shown in Fig. 9(c). The average gray value sum of three intensity patterns is calculated and the total of photons is estimated as 8.77 × 10⁷. Looking up the line plot in Fig. 2(b), the shot noise level is around 3.17%. A stereo display of pattern at the imaging distance of 70mm is also shown in Fig. 9(c). As we can see, there are many burrs in it. The additive noise introduced by dark frame subtraction and alien scattering as well as the quantization error are naturally inherent in the measurements. The missing central intensity problem is not taken into consideration here. Thus, the patterns are mainly corrupted by the shot noise. There are two reasons why the dark illumination is set intentionally. On the one hand, the high-frequency details in diffraction patterns can be kept. On the other hand, it simulates the high shot noise scenario in X-ray experiments that diffraction data at high scattering angles usually has a low signal-to-noise level.

The algorithm parameter settings are kept the same as in the simulation except the support size is 650 × 650 pixels. All three measured noisy patterns are input into multi-image methods while the single-shot methods only utilize the pattern measured at the imaging distance of 70mm.

Figure 9(d) displays the reconstructed results of five algorithms after 5000 iterations. Clearly seen is that SAPR achieves the result closest to the sample shown in Fig. 9(b). The results from APR and NRAPR are acceptable and superior to those from NRHIO and OSS. The experimental results coincide with the conclusions above, considering the dominance of the shot noise. Also, it is shown that the multi-image methods generally outperform the single-shot methods in retrieval fidelity. This is due to more phase information encoded in multiple measurements, thus providing more constraints during the iteration.

In the experiments, E_real is hardly to be depended on, since the true object distribution is almost impossible to acquire in advance. On another hand, E_reciprocal is proved not a reliable indicator of iteration in the subsections above. Therefore, a new metric is proposed to monitor the convergence and the iterative calculation. Inspired by the auto-focusing (AF) process, searching for a global minimum of solution space in phase retrieval parallels locating the focal plane in AF. A well-retrieved image usually has sharp edges and clear detail features, corresponding to the ample high frequency components. Thus, we define a squared-gradient measure

SoG (ψ^{(n)}) = \frac{1}{M N} \sum_{m = 1}^{M} \sum_{n = 1}^{N} [{grad}_{x}^{2} ({| ψ^{(n)} (m, n) |}^{2}) + {grad}_{y}^{2} ({| ψ^{(n)} (m, n) |}^{2})] .

The larger SoG is, the more high-frequency information it contains. The iteration of five algorithms to cope with the experiment data is monitored by SoG, shown in Fig. 10. First of all, the final SoG value of five curves fits well with the observation in Fig. 9(d). Also, the metric ascends gradually as the iteration goes on. Notable is the overshooting peak of SAPR and OSS curves in the first several iterations, which is attributed to the initial guess filtered by the smoothness constraint, shown in the dash box of Fig. 10. But it does not affect the general curve trend. Moreover, the ‘stepping’ phenomenon of SAPR and OSS curves also occurs in E_real curves from Fig. 4, thus validating the effectiveness of SoG. To conclude, the experimental results also demonstrate the superiority of our proposed methods.

Fig. 10 The new defined error metric SoG as function of iteration. The picture in the dash box is the amplitude of $ψ^{(n)}$ at the second iteration of SAPR algorithm.

Download Full Size | PDF

5. Conclusion

In summary, we proposed two noise-robust multiple-intensity phase retrieval algorithms based on the Pauta criterion and the smoothness constraint. Three typical experimental noises are modeled. Then, three old algorithms and two proposed algorithms are investigated under different noises in a broad sense, including the quantization error and the missing central intensity problem. For the shot noise, typical of multiplicative noise, both NRAPR and SAPR have higher retrieval accuracy and convergence speed. For two kinds of additive noise and the quantization error, SAPR has a superior performance at the low noise level and NRAPR is the optimum algorithm on the whole. Confronted with the missing central intensity problem, two proposed algorithms still obtain better reconstructions with less density oscillation. At last, we demonstrate the superiority of proposed algorithms in the experiment and put forward a more reliable and practical metric for experiments, compared with the conventional iteration metrics.

Due to the structure similarity, the modification approach can also be applied to other multiple-intensity phase retrieval methods, like the SBMIR algorithm and the multi-stage algorithm. Recently, these multi-image phase retrieval algorithms have been successfully applied in wavefront sensing [46], numerical aberration correction [47], deformation analysis [48], encryption [49, 50] and pathology diagnosis [51, 52]. The notable improvement of noise robustness demonstrated by numerical analysis and experimental results is quite encouraging, as it raises the prospect of achieving higher resolutions with noisy diffraction data in these applications. It is expected that our work can further expand the application of CDI techniques and iterative phase retrieval methods.

Funding

This work was supported by the National Natural Science Foundation of China (Nos. 61377016, 61575055 and 61575053), the Fundamental Research Funds for the Central Universities (No. HIT.BRETIII.201406), the Program for New Century Excellent Talents in University (No. NCET-12-0148), and the China Postdoctoral Science Foundation (Nos. 2013M540278 and 2015T80340), and the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, China.

Acknowledgments

The authors thank Su Zhang, Cheng Guo and Xue Wang for experiment discussion and the anonymous reviewers for their helpful suggestions and comments.

References and links

1. H. N. Chapman, A. Barty, S. Marchesini, A. Noy, S. P. Hau-Riege, C. Cui, M. R. Howells, R. Rosen, H. He, J. C. H. Spence, U. Weierstall, T. Beetz, C. Jacobsen, and D. Shapiro, “High-resolution ab initio three-dimensional X-ray diffraction microscopy,” J. Opt. Soc. Am. A 23(5), 1179–1200 (2006). [CrossRef] [PubMed]

2. P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer, “High-resolution scanning X-ray diffraction microscopy,” Science 321(5887), 379–382 (2008). [CrossRef] [PubMed]

3. C. T. Putkunz, A. J. D’Alfonso, A. J. Morgan, M. Weyland, C. Dwyer, L. Bourgeois, J. Etheridge, A. Roberts, R. E. Scholten, K. A. Nugent, and L. J. Allen, “Atom-scale ptychographic electron diffractive imaging of boron nitride cones,” Phys. Rev. Lett. 108(7), 073901 (2012). [CrossRef] [PubMed]

4. A. Szameit, Y. Shechtman, E. Osherovich, E. Bullkich, P. Sidorenko, H. Dana, S. Steiner, E. B. Kley, S. Gazit, T. Cohen-Hyams, S. Shoham, M. Zibulevsky, I. Yavneh, Y. C. Eldar, O. Cohen, and M. Segev, “Sparsity-based single-shot subwavelength coherent diffractive imaging,” Nat. Mater. 11(5), 455–459 (2012). [CrossRef] [PubMed]

5. D. F. Gardner, M. Tanksalvala, E. R. Shanblatt, X. Zhang, B. R. Galloway, C. L. Porter, R. Karl Jr, C. Bevis, D. E. Adams, H. C. Kapteyn, M. M. Murnane, and G. F. Mancini, “Subwavelength coherent imaging of periodic samples using a 13.5 nm tabletop high-harmonic light source,” Nat. Photonics 11(4), 259–263 (2017). [CrossRef]

6. J. Miao, P. Ercius, and S. J. L. Billinge, “Atomic electron tomography: 3D structures without crystals,” Science 353(6306), aaf2157 (2016). [CrossRef] [PubMed]

7. P. Thibault and A. Menzel, “Reconstructing state mixtures from diffraction measurements,” Nature 494(7435), 68–71 (2013). [CrossRef] [PubMed]

8. M. M. Seibert, T. Ekeberg, F. R. N. C. Maia, M. Svenda, J. Andreasson, O. Jönsson, D. Odić, 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. Abergel, V. Seltzer, J. M. Claverie, C. Bostedt, J. D. Bozek, S. Boutet, A. A. Miahnahri, 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. Frank, 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. Hömke, C. Reich, D. Pietschner, G. Weidenspointner, L. Strüder, G. Hauser, H. Gorke, J. Ullrich, I. Schlichting, S. Herrmann, G. Schaller, F. Schopper, H. Soltau, K. U. Kühnel, R. Andritschke, C. D. Schröter, 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(7332), 78–81 (2011). [CrossRef] [PubMed]

9. I. Schlichting and J. Miao, “Emerging opportunities in structural biology with X-ray free-electron lasers,” Curr. Opin. Struct. Biol. 22(5), 613–626 (2012). [CrossRef] [PubMed]

10. R. Horstmeyer, J. Chung, X. Ou, G. Zheng, and C. Yang, “Diffraction tomography with Fourier ptychography,” Optica 3(8), 827–835 (2016). [CrossRef]

11. 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(6742), 342–344 (1999). [CrossRef]

12. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Jena) 35, 237–246 (1972).

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

14. J. R. Fienup, “Phase retrieval algorithms: a personal tour [Invited],” Appl. Opt. 52(1), 45–56 (2013). [CrossRef] [PubMed]

15. H. Qi, L. Ruan, M. Shi, W. An, and H. Tan, “Application of multi-phase particle swarm optimization technique to inverse radiation problem,” J. Quant. Spectrosc. Ra. 109(3), 476–493 (2008). [CrossRef]

16. C. Zhang and X. Jian, “Wide-spectrum reconstruction method for a birefringence interference imaging spectrometer,” Opt. Lett. 35(3), 366–368 (2010). [CrossRef] [PubMed]

17. T. Mu, C. Zhang, C. Jia, and W. Ren, “Static hyperspectral imaging polarimeter for full linear Stokes parameters,” Opt. Express 20(16), 18194–18201 (2012). [CrossRef] [PubMed]

18. J. M. Rodenburg and H. M. L. Faulkner, “A phase retrieval algorithm for shifting illumination,” Appl. Phys. Lett. 85(20), 4795–4797 (2004). [CrossRef]

19. A. M. Maiden and J. M. Rodenburg, “An improved ptychographical phase retrieval algorithm for diffractive imaging,” Ultramicroscopy 109(10), 1256–1262 (2009). [CrossRef] [PubMed]

20. G. Zheng, R. Horstmeyer, and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nat. Photonics 7(9), 739–745 (2013). [CrossRef] [PubMed]

21. S. Marchesini, A. Schirotzek, C. Yang, H. Wu, and F. Maia, “Augmented projections for ptychographic imaging,” Inverse Probl. 29(11), 115009 (2013). [CrossRef]

22. T. M. Godden, R. Suman, M. J. Humphry, J. M. Rodenburg, and A. M. Maiden, “Ptychographic microscope for three-dimensional imaging,” Opt. Express 22(10), 12513–12523 (2014). [CrossRef] [PubMed]

23. S. Marchesini, H. Krishnan, B. J. Daurer, D. Shapiro, T. Perciano, J. A. Sethian, and F. R. N. C. Maia, “SHARP: a distributed, GPU-based ptychographic solver,” J. Appl. Cryst. 49(4), 1245–1252 (2016). [CrossRef]

24. A.-L. Robisch, K. Kröger, A. Rack, and T. Salditt, “Near-field ptychography using lateral and longitudinal shifts,” New J. Phys. 17(7), 073033 (2015). [CrossRef]

25. P. Thibault and M. Guizar-Sicairos, “Maximum-likelihood refinement for coherent diffractive imaging,” New J. Phys. 14(6), 063004 (2012). [CrossRef]

26. T. B. Edo, F. Sweeney, C. Lui, and J. M. Rodenburg, “Noise limit on practical electron ptychography,” J. Phys. Conf. Ser. 241(1), 012065 (2010). [CrossRef]

27. L. H. Yeh, J. Dong, J. Zhong, L. Tian, M. Chen, G. Tang, M. Soltanolkotabi, and L. Waller, “Experimental robustness of Fourier ptychography phase retrieval algorithms,” Opt. Express 23(26), 33214–33240 (2015). [CrossRef] [PubMed]

28. C. Zuo, J. Sun, and Q. Chen, “Adaptive step-size strategy for noise-robust Fourier ptychographic microscopy,” Opt. Express 24(18), 20724–20744 (2016). [CrossRef] [PubMed]

29. G. Pedrini, W. Osten, and Y. Zhang, “Wave-front reconstruction from a sequence of interferograms recorded at different planes,” Opt. Lett. 30(8), 833–835 (2005). [CrossRef] [PubMed]

30. J. A. Rodrigo, H. Duadi, T. Alieva, and Z. Zalevsky, “Multi-stage phase retrieval algorithm based upon the gyrator transform,” Opt. Express 18(2), 1510–1520 (2010). [CrossRef] [PubMed]

31. Z. Liu, C. Guo, J. Tan, Q. Wu, L. Pan, and S. Liu, “Iterative phase amplitude retrieval from multiple images in gyrator domains,” J. Opt. 17(2), 025701 (2015). [CrossRef]

32. V. Elser, “Solution of the crystallographic phase problem by iterated projections,” Acta Crystallogr. A 59(3), 201–209 (2003). [CrossRef] [PubMed]

33. D. R. Luke, “Relaxed averaged alternating reflections for diffraction imaging,” Inv. Probl. 21(1), 37–50 (2004). [CrossRef]

34. A. Martin, F. Wang, N. Loh, T. Ekeberg, F. Maia, M. Hantke, G. van der Schot, C. Hampton, R. Sierra, A. Aquila, S. Bajt, M. Barthelmess, C. Bostedt, J. Bozek, N. Coppola, S. Epp, B. Erk, H. Fleckenstein, L. Foucar, M. Frank, H. Graafsma, L. Gumprecht, A. Hartmann, R. Hartmann, G. Hauser, H. Hirsemann, P. Holl, S. Kassemeyer, N. Kimmel, M. Liang, L. Lomb, S. Marchesini, K. Nass, E. Pedersoli, C. Reich, D. Rolles, B. Rudek, A. Rudenko, J. Schulz, R. Shoeman, H. Soltau, D. Starodub, J. Steinbrener, F. Stellato, L. Strüder, J. Ullrich, G. Weidenspointner, T. White, C. Wunderer, A. Barty, I. Schlichting, M. Bogan, and H. Chapman, “Noise-robust coherent diffractive imaging with a single diffraction pattern,” Opt. Express 20(15), 16650–16661 (2012). [CrossRef]

35. 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. Cryst. 46(2), 312–318 (2013). [CrossRef] [PubMed]

36. Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, “Phase retrieval with application to optical imaging,” IEEE Signal Process. Mag. 32(3), 87–109 (2015). [CrossRef]

37. C. Shen, J. Tan, C. Wei, and Z. Liu, “Coherent diffraction imaging by moving a lens,” Opt. Express 24(15), 16520–16529 (2016). [CrossRef] [PubMed]

38. K. A. Nugent, “Coherent methods in the X-ray sciences,” Adv. Phys. 59(1), 1–99 (2010). [CrossRef]

39. Y. Takahashi, A. Suzuki, N. Zettsu, Y. Kohmura, Y. Senba, H. Ohashi, K. Yamauchi, and T. Ishikawa, “Towards high-resolution ptychographic X-ray diffraction microscopy,” Phys. Rev. B 83(21), 214109 (2011). [CrossRef]

40. P. F. Almoro and S. G. Hanson, “Object wave reconstruction by speckle illumination and phase retrieval,” J. Eur. Opt. Soc. 4(4), 09002 (2009). [CrossRef]

41. P. F. Almoro, G. Pedrini, P. N. Gundu, W. Osten, and S. G. Hanson, “Phase microscopy of technical and biological samples through random phase modulation with a diffuser,” Opt. Lett. 35(7), 1028–1030 (2010). [CrossRef] [PubMed]

42. J. Miao, Y. Nishino, Y. Kohmura, B. Johnson, C. Song, S. H. Risbud, and T. Ishikawa, “Quantitative image reconstruction of GaN quantum dots from oversampled diffraction intensities alone,” Phys. Rev. Lett. 95(8), 085503 (2005). [CrossRef] [PubMed]

43. http://www.imageprocessingplace.com/root_files_V3/image_databases.htm

44. G. Williams, M. Pfeifer, I. Vartanyants, and I. Robinson, “Effectiveness of iterative algorithms in recovering phase in the presence of noise,” Acta Crystallogr. A 63(1), 36–42 (2007). [CrossRef] [PubMed]

45. J. Miao and D. Sayre, “On possible extensions of X-ray crystallography through diffraction-pattern oversampling,” Acta Crystallogr. A 56(6), 596–605 (2000). [CrossRef] [PubMed]

46. P. F. Almoro and S. G. Hanson, “Wavefront sensing using speckles with fringe compensation,” Opt. Express 16(11), 7608–7618 (2008). [CrossRef] [PubMed]

47. P. F. Almoro, P. N. Gundu, and S. G. Hanson, “Numerical correction of aberrations via phase retrieval with speckle illumination,” Opt. Lett. 34(4), 521–523 (2009). [CrossRef] [PubMed]

48. A. Anand, V. K. Chhaniwal, P. Almoro, G. Pedrini, and W. Osten, “Shape and deformation measurements of 3D objects using volume speckle field and phase retrieval,” Opt. Lett. 34(10), 1522–1524 (2009). [CrossRef] [PubMed]

49. G. Cheng, C. Wei, J. Tan, K. Chen, S. Liu, Q. Wu, and Z. Liu, “A review of iterative phase retrieval for measurement and encryption,” Opt. Lasers Eng. 89, 2–12 (2016).

50. Z. Liu, C. Shen, J. Tan, and S. Liu, “A recovery method of double random phase encoding system with a parallel phase retrieval,” IEEE Photonics J. 8(1), 7801807 (2016). [CrossRef]

51. A. Greenbaum, Y. Zhang, A. Feizi, P. L. Chung, W. Luo, S. R. Kandukuri, and A. Ozcan, “Wide-field computational imaging of pathology slides using lens-free on-chip microscopy,” Sci. Transl. Med. 6(267), 267ra175 (2014). [CrossRef] [PubMed]

52. W. Luo, A. Greenbaum, Y. Zhang, and A. Ozcan, “Synthetic aperture-based on-chip microscopy,” Light Sci. Appl. 4(3), e261 (2015). [CrossRef]

Two noise-robust axial scanning multi-image phase retrieval algorithms based on Pauta criterion and smoothness constraint

Abstract

1. Introduction

2. Algorithm analysis

3. Noise model

4. Results and discussion

4.1 Photon number

4.2 Dark frame subtraction

4.3 Alien scattering

4.4 Quantization error

4.5 Center missing

4.6 Experimental demonstration

5. Conclusion

Funding

Acknowledgments

References and links

Cited By

Figures (10)

Tables (2)

Equations (6)

Optics Express