## Abstract

We propose a novel and effective method to quantitatively recover a complex-valued object from diffraction intensity maps recorded in the fractional Fourier domain. A wavefront modulation is introduced in the wave path, and several diffraction intensity maps are recorded through variable function orders in the fractional Fourier transform. A new phase retrieval algorithm is then proposed, and advantages of the proposed algorithm are also discussed. A proof-of-principle study is presented to show the feasibility and effectiveness of the proposed method.

© 2010 OSA

## 1. Introduction

Phase retrieval technique plays an important role in the structural determination of micro or nano-objects [1–3]. Because only intensity distribution can be recorded by a charge-coupled device (CCD) camera, the phase information is lost. Since holography was introduced by Gabor [3], it has been widely applied to various fields [4–6], such as biological microscopy [4] and metrology [5,6]. Two main types of holography, i.e., off-axis digital holography [4] and phase-shifting digital holography [7], have been developed. Since complex amplitude can be extracted by a numerical reconstruction algorithm in digital holography, phase information is directly determined by using an arc-tangent operation. However, the interferometric method requires a relatively complex optical system, such as a reference wave, temporal coherence and a stable mechanical setup. Recently, non-interferometric method [1] is considered as a promising alternative for phase retrieval, and has attracted much attention in various fields, such as material [1,2].

Hybrid input-output method which was partially derived based on Gerchberg-Saxton algorithm [8] was proposed by Fienup [9], and it has been commonly used in the non-interferometric phase retrieval. When the test object is real or non-negative, it is easy to apply certain support constraints to extract the phase information from the recorded diffraction intensity. However, support constraints in the conventional hybrid input-output method [9] might be inapplicable for the study of complex-valued objects, and stagnation problem is still a great challenge. It is demonstrated that a stronger support constraint is usually required for a complex-valued object than that for a real or non-negative object. Miao et al. [10] proposed an oversampling phasing method to extract complex-valued objects. Liu et al. [11] proposed a phase contrast imaging technique based on a combination of structured illumination and an optimized hybrid input-output algorithm for the extraction of complex-valued objects. In this paper, we propose a novel method with variable function orders in the fractional Fourier transform (FRFT) for the extraction of a complex-valued object. A wavefront modulation is introduced in the wave path, and several diffraction intensity maps are recorded. A new phase retrieval algorithm is then proposed to recover the complex-valued object.

## 2. Theoretical analysis

Figure 1 shows a schematic experimental arrangement. A collimated plane wave is used to illuminate a test complex-valued object, and a diffraction intensity distribution is recorded by a CCD camera. A series of intensity distributions can be recorded, when the function order in the FRFT is modified correspondingly. Wave propagation between the object plane and the phase-only mask plane can be described by the diffraction principle.

*λ*is the wavelength of a light source, wave number $k=2\pi /\lambda ,$

*ρ*is the distance from a point in the object plane to a point in the phase-only mask plane$\left[\rho =\sqrt{{(x-\xi )}^{2}+{(y-\eta )}^{2}+{z}^{2}}\right],$ and

*z*denotes the axial distance between the object plane and the phase-only mask plane.

Equation (1) can be simplified with a paraxial or small-angle approximation, and can also be implemented by a convolution method. In this study, angular spectrum algorithm [5,12] is used, and the complex amplitude just before the phase-only mask plane is described by

The object wave *F(x,y)* is multiplied by a known phase-only mask *M(x,y)*, and then propagates through a fractional Fourier domain. In practice, a spatial light modulator can be used to display the phase-only mask *M(x,y)*. For the sake of brevity, an one-dimensional FRFT [13–16] is described, and the FRFT with an order *a* can be defined as

*m*denotes an integer, and $R=\sqrt{1-j\text{\hspace{0.05em} \hspace{0.05em}}\mathrm{cot}\left(a\pi /2\right)}.$ Similarly, the 2D FRFT can be defined, and an intensity map is expressed aswhere for simplicity the value of

*b*is equal to

*a*in this study. The function order can be modified through the adjustment of the distance and the focal length in the FRFT domain [13–16], thus a series of diffraction intensity maps are recorded correspondingly.

Here, we propose a new iterative algorithm using three recorded diffraction intensity maps $\left[{I}^{1}(\mu ,\nu ),\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{I}^{2}(\mu ,\nu )\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}\mathrm{and}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{I}^{3}(\mu ,\nu )\right]$ to recover the original complex-valued object, and a flow chart for the proposed method is shown in Fig. 2 . The proposed phase retrieval algorithm proceeds as follows:

- (1) Assume an initial random or constant distribution ${F}_{n}(x,y)$ (iteration number n = 1, 2, 3…) just before the phase-only mask $M(x,y);$
- (2) Multiply by the phase-only mask$M(x,y):$ ${F}_{n}^{\text{'}}(x,y)={F}_{n}(x,y)\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}M(x,y);$

- (3) Propagate to CCD plane with function order of ${a}_{1}:$ ${O}_{n}^{1}(\mu ,v)={\mathrm{FRFT}}_{{a}_{1},{a}_{1}}\text{\hspace{0.05em}}\left[{F}_{n}^{\text{'}}(x,y)\right]\text{\hspace{0.05em} \hspace{0.05em}};$
- (4) Apply a support constraint in the CCD plane with the square root of the recorded intensity distribution: $\overline{{O}_{n}^{1}(\mu ,\nu )}={\left[{I}^{1}(\mu ,\nu )\right]}^{1/2}\left[{O}_{n}^{1}(\mu ,\nu )/\left|{O}_{n}^{1}(\mu ,\nu )\right|\right]\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em}};$
- (5) Numerically propagate back to the phase mask plane and multiply by${M}^{*}(x,y):$ $\overline{{F}_{n}^{1}(x,y)}={\mathrm{FRFT}}_{-{a}_{1},-{a}_{1}}\left[\overline{{O}_{n}^{1}(\mu ,\nu )}\right]\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}{M}^{*}(x,y),$ and asterisk denotes complex conjugate;
- (6) Using $\overline{{F}_{n}^{1}(x,y)},$ repeat steps 2-5 with corresponding parameter modifications (i.e., function order and recorded intensity map), and obtain $\overline{{F}_{n}^{2}(x,y)};$
- (7) Using $\overline{{F}_{n}^{2}(x,y)},$ further repeat steps 2-5 with corresponding parameter modifications, and obtain $\overline{{F}_{n}^{3}(x,y)};$
- (8) Calculate the measurement error by using the following formulae:$${{\displaystyle \sum \text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}\left[\text{\hspace{0.05em} \hspace{0.05em}}\left|\text{\hspace{0.05em}}\overline{{F}_{n}^{3}(x,y)}\text{\hspace{0.05em}}\right|-\left|{F}_{n}(x,y)\right|\right]}}^{2}\le \sigma \text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}(if\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}n=1);$$$${{\displaystyle \sum \text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}\left[\text{\hspace{0.05em} \hspace{0.05em}}\left|\text{\hspace{0.05em}}\overline{{F}_{n}^{3}(x,y)}\text{\hspace{0.05em}}\right|-\text{\hspace{0.05em}}\left|\text{\hspace{0.05em}}\overline{{F}_{n-1}^{3}(x,y)}\text{\hspace{0.05em}}\right|\right]}}^{2}\le \sigma \text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}(\mathit{\text{if}}\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}n>1).$$

If the measurement error is not larger than a preset threshold (*σ*), the iteration operation stops. Otherwise, the complex amplitude $\overline{{F}_{n}^{3}(x,y)}$ obtained in the step 7 is considered as the object wave just before the phase-only mask plane for the next iteration (n = n + 1).

- (9) When the condition is met, the complex-valued object can be recovered by using the angular spectrum algorithm with an inverse distance (-
*z*).

To evaluate the accuracy of the recovered object, sum-squared error (SSE) is calculated by

## 3. Results and discussion

We carry out a numerical experiment as shown in Fig. 1 to demonstrate the feasibility and effectiveness of the proposed method. The pixel size of CCD camera is $4.65\text{\hspace{0.05em} \hspace{0.05em} \hspace{0.05em}}\mu m,$ and the pixel number is$512\times 512$. The light source wavelength is 632.8 nm, and the distance between the object plane and phase-only mask plane is 20 mm. The phase mask is a map randomly distributed in a range of 0 to$2\pi .$ According to the definition of FRFT [13–16], function orders of 0.25, 0.50 and 0.75 are generated in the three recordings, respectively. Original object amplitude (Barbara) and phase (Phantom) are shown in Figs. 3(a) and 3(b), and a typical one of recorded diffraction intensity maps is shown in Fig. 3(c). With the proposed phase retrieval algorithm, object amplitude and phase are recovered as shown in Figs. 3(d) and 3(e). The phase-contrast map is directly extracted by an arc-tangent operation from the recovered complex-valued object. The number of iterations is 63 to ensure that measurement error is not larger than a preset threshold of $\sigma =\mathrm{0.001.}$ As the threshold increases, the number of iterations will reduce correspondingly and the quality of reconstructed images might degrade. To present the quality of the recovered phase map, a comparison along the dashed line [as indicated in Fig. 3(b)] between Figs. 3(b) and 3(e) is shown in Fig. 3(f). It can be seen that a high-quality phase map is recovered, and the only difference from the original phase map is just a constant offset. Figure 4(a) shows a relationship between the number of iterations and SSE. It can be seen that the SSE decreases rapidly just after several iterations.

The study further demonstrates that the number of iterations decreases considerably with more recordings, but more experimental effort should be made. For instance, with 4 recordings, the iteration number decreases to 27 since more information related to the object is available. Figure 4(b) shows the relationship between the number of iterations and SSE using four recorded intensity maps. It is noteworthy that besides an additional function order (0.95), other experimental parameters are the same as those in Fig. 3. The proposed method with only two diffraction intensity recordings is feasible in theory, but the quality of the reconstructed images may not be good and some distortions could make the iterative operation stagnant.

Since the recordings are usually disturbed by noise, the recorded diffraction intensity maps with different levels of noise are also studied. Figures 5(a) and 5(b) show the relationships between the number of iterations and SSE using three and four noisy intensity maps, respectively. Random noise with a signal-to-noise ratio of 20 is added to all intensity maps, and the number of iterations is set as 300. The fast convergence rate is also achieved, and the proposed method still can effectively extract the information about the complex-valued object but probably with a relatively low quality of the extracted phase-contrast map and an increase of the number of iterations. In this case, an average algorithm can be applied to a series of recovered phase maps using different initial estimates.

In the conventional phase-retrieval algorithm, the twin image may exist which will highly affect the quality of the recovered images. In the proposed method, a phase mask with randomly-distributed values is placed in the wave path, which results in more uncorrelated diffraction distributions. Therefore, the twin image problem can be addressed, and the stagnation problem is effectively mitigated. Note that the phase mask is not modified during the recordings, and no constraint is applied to the phase mask during phase retrieval. It is illustrated that the proposed method can effectively and accurately extract the information about the complex-valued object. However, the convergence rate and the algorithm stagnation could be related to how large the difference among the variable function orders is.

## 4. Conclusions

In this paper, we have proposed a novel and effective method to quantitatively recover a complex-valued object from the diffraction intensity maps. A series of diffraction intensity maps are recorded by using variable function orders in the FRFT, and a new phase retrieval algorithm is then proposed. It is illustrated that the recovered phase and amplitude images are of high quality. In addition, a fast convergence rate is achieved, and no support constraint is required in the phase mask. The proposed method can also be applied in the extended FRFT domain [17] and find some other applications, such as in the virtual optics.

## Acknowledgements

This work was supported by the Singapore Temasek Defence Systems Institute under grant TDSI/09-001/1A, and the authors are grateful to the reviewers for their valuable comments.

## References and links

**1. **J. M. Zuo, I. Vartanyants, M. Gao, R. Zhang, and L. A. Nagahara, “Atomic resolution imaging of a carbon nanotube from diffraction intensities,” Science **300**(5624), 1419–1421 (2003). [CrossRef]

**2. **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**(7098), 63–66 (2006). [CrossRef]

**3. **D. Gabor, “A new microscopic principle,” Nature **161**(4098), 777–778 (1948). [CrossRef]

**4. **C. J. Mann, L. Yu, C. Lo, and M. K. Kim, “High-resolution quantitative phase-contrast microscopy by digital holography,” Opt. Express 13, 8693–8698 (2005), http://www.opticsinfobase.org/oe/abstract.cfm? URI = oe-13–22–8693.

**5. **W. Chen, C. Quan, and C. J. Tay, “Extended depth of focus in a particle field measurement using a single-shot digital hologram,” Appl. Phys. Lett. **95**(20), 201103 (2009). [CrossRef]

**6. **W. Chen, C. Quan, and C. J. Tay, “Measurement of curvature and twist of a deformed object using digital holography,” Appl. Opt. **47**(15), 2874–2881 (2008). [CrossRef]

**7. **I. Yamaguchi and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. **22**(16), 1268–1270 (1997). [CrossRef]

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

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

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

**11. **Y. J. Liu, B. Chen, E. R. Li, J. Y. Wang, A. Marcelli, S. W. Wilkins, H. Ming, Y. C. Tian, K. A. Nugent, P. P. Zhu, and Z. Y. Wu, “Phase retrieval in X-ray imaging based on using structured illumination,” Phys. Rev. A **78**(2), 023817 (2008). [CrossRef]

**12. **J. W. Goodman, *Introduction to Fourier Optics*, 2nd ed. (McGraw-Hill, 1996).

**13. **A. W. Lohmann, “Image rotation, Wigner rotation, and the fractional Fourier transform,” J. Opt. Soc. Am. A **10**(10), 2181–2186 (1993). [CrossRef]

**14. **H. M. Ozaktas, Z. Zalevsky, and M. A. Kutay, *The Fractional Fourier Transform with Applications in Optics and Signal Processing* (Wiley, 2001).

**15. **M. G. Ertosun, H. Atlı, H. M. Ozaktas, and B. Barshan, “Complex signal recovery from two fractional Fourier transform intensities: order and noise dependence,” Opt. Commun. **244**(1-6), 61–70 (2005). [CrossRef]

**16. **M. G. Ertosun, H. Atlı, H. M. Ozaktas, and B. Barshan, “Complex signal recovery from multiple fractional Fourier-transform intensities,” Appl. Opt. **44**(23), 4902–4908 (2005). [CrossRef]

**17. **J. Hua, L. Liu, and G. Li, “Extended fractional Fourier transforms,” J. Opt. Soc. Am. A **14**(12), 3316–3322 (1997). [CrossRef]