Multi-stage phase retrieval algorithm based upon the gyrator transform

The gyrator transform is a useful tool for optical information processing applications. In this work we propose a multi-stage phase retrieval approach based on this operation as well as on the well-known Gerchberg-Saxton algorithm. It results in an iterative algorithm able to retrieve the phase information using several measurements of the gyrator transform power spectrum. The viability and performance of the proposed algorithm is demonstrated by means of several numerical simulations and experimental results. © 2009 Optical Society of America OCIS codes: (070.2590) Fourier transforms; (120.4820) Optical systems; (200.4740) Optical processing; (140.3300) Laser beam shaping References 1. D. Gabor, “A new microscopic principle,” Nature 161, 777–778 (1948). 2. E. N. Leith and J. Upatnieks, “Reconstructed wavefronts and communication theory,” J. Opt. Soc. Am. 52, 1123– 1128 (1962). 3. E. N. Leith and J. Upatnieks, “Wavefront reconstruction with continuous-tone objects,” J. Opt. Soc. Am. 53, 1377–1381 (1963). 4. E. N. Leith and J. Upatnieks, “Wavefront reconstruction with diffused illumination and three-dimensional objects,” J. Opt. Soc. Am. 54, 1295–1301 (1964). 5. J. P. Guigay, “Fourier-transrorm analysis of Fresnel diffraction patterns and in-line holograms,” Optik 49, 121– 125 (1977). 6. M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. 73, 1434–1441 (1983). 7. T. E. Gureyev and K. A. Nugent, “Rapid quantitative phase imaging using the transport of intensity equation,” Opt. Commun. 133, 339 – 346 (1997). 8. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972). 9. M. Nieto-Vesperinas, R. Navarro, and F. J. Fuentes, “Performance of a simulated-annealing algorithm for phase retrieval,” J. Opt. Soc. Am. A 5, 30–38 (1988). 10. J. H. Seldin and J. R. Fienup, “Iterative blind deconvolution algorithm applied to phase retrieval,” J. Opt. Soc. Am. A 7, 428–433 (1990). 11. J. R. Fienup, “Phase-retrieval algorithms for a complicated optical system,” Appl. Opt. 32, 1737–1746 (1993). 12. D. L. Misell, “A method for the solution of the phase problem in electron microscopy,” J. Phys. D: Applied Physics 6, L6–L9 (1973). 13. D. C. Redding, S. Basinger, A. Lowman, A. Kissil, P. Bely, R. Burg, R. Lyon, G. Mosier, M. Femiano, M. Wilson, R. G. Schunk, L. Craig, D. Jacobson, J. Rakoczy, and J. Hadaway, “Wavefront Sensing and Control for a NextGeneration Space Telescope,” Proc. SPIE 3356, 758 (1998). #118097 $15.00 USD Received 2 Oct 2009; revised 16 Nov 2009; accepted 23 Nov 2009; published 12 Jan 2010 (C) 2010 OSA 18 January 2010 / Vol. 18, No. 2 / OPTICS EXPRESS 1510 14. D. S. Acton, P. D. Atcheson, M. Cermak, L. K. Kingsbury, F. Shi, and D. C. Redding, “James Webb Space Telescope wavefront sensing and control algorithms,” Proc. SPIE 5487, 887 (2004). 15. G. R. Brady and J. R. Fienup, “Nonlinear optimization algorithm for retrieving the full complex pupil function,” Opt. Express 14, 474–486 (2006). 16. B. H. Dean, D. L. Aronstein, J. S. Smith, R. Shiri, and D. S. Acton, “Phase retrieval algorithm for JWST Flight and Testbed Telescope,” Proc. SPIE 6265, 626511 (2006). 17. Z. Zalevsky, D. Mendlovic, and R. G. Dorsch, “Gerchberg-Saxton algorithm applied in the fractional Fourier or the Fresnel domain,” Opt. Lett. 21, 842 (1996). 18. J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Experimental implementation of the gyrator transform,” J. Opt. Soc. Am. A 24, 3135–3139 (2007). 19. J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Optical system design for orthosymplectic transformations in phase space,” J. Opt. Soc. Am. A 23, 2494–2500 (2006).


Introduction
Phase retrieval plays an important role in science, engineering, and in particular in optical imaging and information processing.In general, its experimental application demands a flexible and robust setup which is an important issue.There exist several approaches for phase retrieval used for example in the estimation of the depth dimension of objects, characterization of optical elements, etc.Since holography permits the recording of the complex field amplitude, some of these phase retrieval approaches exploit interferometric techniques [1,2,3,4], which requires the use of a reference signal.Alternatively, phase retrieval approaches based only on the measurement of intensity distributions during signal propagation are applied.
These methods are either iterative or non-iterative involving digital post processing.In particular, the non-iterative approaches are able to retrieve the phase information by means of transforming Guigay equations [5,6], by solving the so-called transport of intensity equation [7], etc.Nevertheless, they are often mathematically complex.Other option is to use iterative techniques based on the well-known Gerchberg-Saxton (GS) algorithm [8].The GS algorithm is widely used in different optical applications and permits to retrieve the input phase from the known amplitude distributions corresponding to input and output planes of the optical setup.Its main problem for phase retrieval is the phase ambiguity resulting with converging to local minima solutions, which yields an incorrect phase distribution.There exist several works trying to solve this problem by using statistical methods [9], but they demand a high computational complexity.Blind de-convolution methods [10] speed up the convergence rate but do not solve the local minima problem either.
Other iterative approaches have been also proposed, for instance the gradient search method [11].However, unlike the GS approach, such techniques are only suitable for special conditions.Alternatively, there exist methods in microscopy using plurality of de-focused images [12,13].These methods use de-focused images as well as a focused image to reconstruct the phase distribution of one of the planes, by using the known relations between the images.Later on, other works were introduced combining this idea with an iterative GS algorithm [14,15,16].This approach was found to be faster than the iterative GS technique and it also solved the phase ambiguity issue due to the usage of multiple constraint planes.Originally GS technique involved a Fourier or Fresnel transform relation between the input and output planes.Nevertheless, later on other transformations relating such planes were suggested [17] as well.
In this work we apply a GS-based approach involving several intensity distributions working as constraint planes.In particular, we consider the gyrator transform (GT) which is mathematically defined as a linear canonical integral transform that produces the twisted rotation in position−spatial frequency planes of phase space [18].Thus, the GT operation of a twodimensional function f i (r i ) (complex field amplitude) with parameter α, known as transforma-tion angle, is given by where r i,o = (x i,o , y i,o ) are the input and output coordinates, respectively.This transformation can be implemented optically by using cylindrical lenses which are properly rotated to change the transformation angle α [19].The experimental implementation of the GT has been demonstrated in Ref. [18].
In our approach at each iteration loop an iterative calculation based on the GT operation for a certain set of angles α is performed.After calculating all the considered angles, an average estimation is obtained yielding a retrieved phase to be used as input for the next iteration loop.The validity of the method is demonstrated by numerical simulations as well as by experimental results.
We remark that GT operation for different angles α can be performed by only the proper rotation of cylindrical lenses, where the distance between them and input-output planes are fixed.This fact avoid several problems arising from misalignment which is an important issue in optical systems that require axial movements.Therefore the proposed approach results attractive for experimental phase retrieval applications.

Multi-stage phase retrieval algorithm
The proposed phase retrieval algorithm is based on the measurements of the intensity distributions associated with GT power spectrum |F α (r o )| 2 for various angles α and Furthermore we will use the following notations for the input signal's amplitude ) and the amplitude of its GT operation , where j = 1, 2, ..., M. Let us now explain the flow chart of the proposed approach as it is depicted in Fig. 1.An initial random phase φ 1 (x, y) ∈ (0, 2π) is chosen.The retrieved phase distribution φ k (x, y) at iteration number k is obtained as follows: 1.The phase estimation from the previous iteration step is multiplied by the input amplitude distribution E in (x, y), resulting with E in (x, y) exp [iφ k (x, y)].

Perform the following calculation M times:
(a) Gyrator transform of E in (x, y) exp [iφ k (x, y)] with the transformation angle α j resulting in a function denoted by: with the transformation angle −α j (which equals to the inverse GT) resulting with A j k (x, y) exp iφ j k (x, y) .
(d) Use the retrieved phase φ j k (x, y) for the next iteration corresponding to angle α j+1 .
The previous calculation leads to a new phase estimation to be used at the next loop: φ k+1 (x, y) = φ M k (x, y).The phase retrieval algorithm ends after N loops yielding a retrieved phase φ N (x, y).Therefore, this approach requires N×M iterations, where M is the number of transformation angles α or considered constrain planes.
In order to demonstrate the phase retrieval capabilities of the proposed approach, we preformed numerical simulations following the parameters used further in the experimental GT Replace amplitude Random phase: Fig. 1.Flow chart corresponding to the proposed algorithm.
setup.In all simulations we have chosen the parameters to match the values of the experimental realization for the GT, which is described by: where λ is the wavelength, z corresponds to fixed axial interval between the lenses and s 2 = 2λ z is a constant coordinate normalization [18].In our experimental setup λ = 532 nm and z = 0.94 cm.
As input signals we consider the Hermite−Gaussian (HG) modes and Laguerre−Gaussian (LG) modes where: r 2 = x 2 + y 2 , H m is the Hermite polynomial, w = √ 2λ z is the beam waist, and L l p is the Laguerre polynomial with radial index p and azimuthal index l.In particular we study the HG 4,3 and LG + 3,1 mode, see Fig. 2. We underline that HG 4,3 is transformed into LG + 3,1 mode under the GT for α = (2q + 1)π/4 and vice versa, where q is an integer, p = min(m, n) and l = |m − n|.
Now we aim to demonstrate that using several (rather than a single) constraint planes for different transformation angles of the GT, produce improved phase retrieval.In particular we study the phase-retrieval performance corresponding to HG 4,3 mode by using one, two and three constraint planes (M = 1, 2, 3...) in which the amplitude distribution is known.The retrieved phase distributions after 30 iterations for HG 4,3 mode are displayed in Fig. 3 The relative root mean square (RMS) error, also known as Abserr, of the retrieved phase at each iteration is given by RMS = ∑ E in e iφ − E in e iφ rt2 1 2 × ∑ E in e iφ 2 − 1 where φ (x, y) is the phase to retrieve and φ rt (x, y) is the estimated phase distribution at each iteration.
The evolution of RMS error for the case of HG 4,3 mode is displayed in Fig. 3(b), respectively calculated for one, two and three constraint planes.We conclude that increasing the number of constraint planes from one to two and three, speeds up the converging of the RMS error as one appreciates in Fig. 3(b).This fact has also been proven for other values of the transformation angles yielding a retrieved phase with a RMS error value of 0.01 or less.
In order to demonstrate the RMS error evolution as a function of the number of constraint planes with low number of iterations, we now consider LG + 3,1 mode as input signal.Specifically, we use only 36 iterations and constraint planes corresponding to the GT operation at angles error value of 0.006 after 36 iterations.In conclusion, when a low number of iterations is utilized more constraint planes are needed in order to obtain an accurate retrieved phase (low RMS value).Meanwhile, as it was previously demonstrated, the RMS error can be reduced to 0.01 or less for the case of only two constraint planes by increasing the number of iterations.To simulate the phase retrieval from measurements of noisy intensity distributions, we consider a normal random additive noise.In particular, we produce this random noise yielding a signal-to-noise ratio (SNR) of 20 dB which is added to each power spectrum or constraint plane.We have chosen this noise and SNR value because it is similar to the one observed in our experimental data.For instance, the constraint plane α = 80 o displayed in Fig. 5(a) is now disturbed by this additive noise as it is shown in Fig. 5(b).Here we perform the same calculation as in the noiseless case: five constraint planes corresponding to angles α = 45 o , 55 o , 65 o , 70 o and 80 o .As it is demonstrated in Fig. 6 the phase distribution is successfully retrieved even in this noisy case by using at least three constraint planes and after 36 iterations (RMS of 0.02).Therefore the proposed algorithm is comparatively robust for noise, SNR of 20 dB in our case.We have also studied the case of high-moderate noisy images.For example in the Table 1 the RMS error values obtained in the case of SNR of 15 dB are also displayed.In contrast to the noiseless case (see RMS error value, as it was expected.We underline that the RMS error is again minimized by increasing the number of iterations as well as constraint planes utilized until the convergence occurs.

Experimental results
The GT setup is constructed by using three generalized lenses (with L 3 = L 1 ) and two fixed intervals z, see Fig. 7(a).Each generalized lens is an assembled set of two thin convergent cylindrical lenses, Fig. 7(b).The rotation angles of the cylindrical lenses are φ 1 = −ϕ 1,2 and φ 2 = −(φ 1 + π/2), where the angles ϕ 1 and ϕ 2 correspond to generalized lens L 1 and L 2 as follows: sin 2ϕ 1 = cot (α/2) /2 and sin 2ϕ 2 = sin α.Notice that we have chosen a coordinate normalization of s 2 = 2λ z instead of the one reported in Ref. [18] (s 2 = λ z).Since this application only demands the acquisition of the intensity distribution then the third generalized lens is not required.Therefore the experimental setup consists of two generalized lenses (four cylindrical lenses), see Fig.We underline that the GT optical setup proposed for the measurements of the required GT power spectra is a flexible and robust system that does not require axial movements.The variation of the GT parameter is achieved by only the proper rotation of cylindrical lenses.These facts make attractive the application of the discussed approach for phase retrieval of optical beams.

Fig. 4 .
Fig. 4. Retrieved phase distribution after N × M = 23 iterations for the case of LG + 3,1 mode, top panel.The phase retrieval is calculated by using constraint planes involving angles: α = 45 o , 55 o , 65 o , 70 o and 80 o .The evolution of RMS error is also displayed as a function of the number of constraint planes as well as for the iterations utilized.

Fig. 5 .(Fig. 6 .
Fig. 5. (a) Intensity distribution (constraint plane) corresponding to the LG + 3,1 mode transformed under GT operation with angle α = 80 o .This image as well as the rest of constraint planes are disturbed by a random additive noise with a SNR value of 20 dB and 15 dB, as displayed in (b) and (c) respectively.It permits to study the proposed phase retrieval algorithm from noisy images.To help visualization these images have been saturated.

Fig. 7 .
Fig. 7. (a) GT optical setup.Each generalized lens L 1,2 is an assembled set of two cylindrical lenses (b).The transformation is only reached by the proper rotation (with angles φ 1 and φ 2 ) of the lenses.

Fig. 8 .
Fig. 8. (a) Experimental GT setup.The intensity distribution of the output signal is registered by a CCD camera with pixel size of 4.6 μm.(b) Generalized lens constructed by assembling two plano-convex cylindrical lenses, which are rotated at angle φ 1 and φ 2 respectively.

8
(a).We remark that each generalized lens is an assembled set of two identical plano-convex cylindrical lenses which are fabricated from N-BK7 #118097 -$15.00USD Received 2 Oct 2009; revised 16 Nov 2009; accepted 23 Nov 2009; published 12 Jan 2010 (C) 2010 OSA algorithm have been applied to the experimental data.The results are in good agreement with the theoretical predictions demonstrating the practical feasibility of the proposed method.

Table 1 .
RMS error values after N ×M = 36 iterations.M is the number of constraint planes.