## Abstract

The standard Gerchberg–Saxton (GS) algorithm is normally used to find the phase (measured on two different parallel planes) of a propagating optical field (usually far-field propagation), given that the irradiance information on those planes is known. This is mostly used to calculate the modulation function of a phase mask so that when illuminated by a plane wave, it produces a known far-field irradiance distribution, or the equivalent, to calculate the phase mask to be used in a Fourier optical system so the desired pattern is obtained on the image plane. There are some extensions of the GS algorithm that can be used when the transformations that describe the optical setup are non-unitary, for example the Yang-Gu algorithm, but these are usually demonstrated using nonunitary translational-invariant optical systems. In this work a practical approach to use the GS algorithm is presented, where raytracing together with the Huygens-Fresnel principle are used to obtain the transformations that describe the optical system, so the calculation can be made when the field is propagated through a translational-variant optical system (TVOS) of arbitrary complexity. Some numerical results are shown for a system where a microscope objective composed by 5 lenses is used.

© 2013 Optical Society of America

## 1. Introduction

The Gerchberg-Saxton (GS) algorithm [1] was proposed in 1972 by R. W. Gerchberg and W. O. Saxton to solve the problem of phase retrieval of a field at two different planes, when at those, only the field amplitudes are known and given that the fields are related by a Fourier transform. Besides its original purpose, it has also been extensively used to calculate phase-only diffractive optical elements (DOE) to be used as beam-shapers [2, 3] or in imaging applications [4], where the optical field’s propagation from the DOE to the image plane and vice versa are known [5] and Fourier like, such as systems involving far-field propagation or in a well-corrected Fourier system (Fig. 1).

There are also some applications related to beam-shaping and imaging, where it is of interest to apply some constrains to the obtained DOE, such as discretizing its phase to simplify the fabrication. Here, some modifications of the standard GS algorithm have been proposed and analyzed [6], but the optical systems involved are usually simulated by a Fourier transform too. Besides, there are some extensions of the GS algorithm that can be used when the transformations that describe the optical setup are non-unitary (and even translational-variant), for example the Yang-Gu algorithm [7], but these are usually demonstrated using simple non-unitary and translational-invariant optical systems such a (paraxial) Fourier system that includes an aperture, and where a method to characterize more complex systems is not addressed. In addition, some other modifications of the GS algorithm where the object (phase mask) and image plane are not related by a Fourier transform but by a Fresnel or a Fractional Fourier transform have been proposed [8], but these algorithms have been developed for and applied to simple systems involving ideal lenses and free space propagation.

A first approximation to handle non-ideal systems, is to regard them as paraxial (translational-invariant) optical systems with aberrations. Here the GS algorithm still works, provided that the adequate compensation phase is added to the calculated phase mask. However, in applications where the optical setup is a translational-variant and non paraxial Fourier-like system that cannot be simulated by ideal lenses and free space propagation, the GS algorithm or any of its derivatives cannot be used, because they do not addresses the problem of how to generate the optical transformations that characterize the system adequately.

What it is proposed in this work is to use a mixed approach, where raytracing and the Huygens-Fresnel principle [9] are used to characterize an arbitrarily complex optical setup, making it possible to obtain the adequate direct and inverse propagation functions to be used with the GS algorithm to calculate the phase mask needed to produce an arbitrary irradiance distribution in the image plane of the setup. As an example, the method is tested using a highly translational-variant multi-lens Fourier-like system.

## 2. Propagation of a field through a TVOS

The first step in determining the propagation of an optical field through a TVOS (translational-variant optical system), is to find a discrete description of it. On a computer this is done when the optical field is represented as a matrix of complex numbers.

The next step is to replace each pixel with a point source, and to obtain by raytracing through the optical system, the (geometric) phase of its propagated field. This is done by assuming the phase difference between a point source on the source-plane and a point (*x*, *y*) on the destination-plane, to be equal to:

*d*is the optical path of the ray joining the point source and the (

*x*,

*y*) point, and

*λ*is its wavelength. One approximation done in this step, is that the amplitude of the field on the observation-plane is constant and proportional to the amplitude of the corresponding source. After the fields due to every source are found, the superposition principle is used and the resulting field is obtained as the sum of them. This process can be performed to propagate the field in any direction through the optical system, so it is possible to propagate the field from the object plane to the image plane, and vice versa.

Mathematically this can be expressed as:

*(*

_{im}*x*,

*y*) is the total optical field at the destination-plane,

*A*and

_{lm}*ϕ*are the amplitude and phase of field on the source-plane at (

_{lm}*x*,

_{l}*y*), and

_{m}*θ*(

_{lm}*x*,

*y*) is the phase of the optical field on the destination-plane due to a point source located at (

*x*,

_{l}*y*) on the source-plane. When propagating forward in Fig. 1, the source-plane is the phase mask and the destination-plane is the image plane. When propagating backward, the source-plane is the image plane, and the destination-plane is the phase mask.

_{m}## 3. The implemented GS algorithm

A block diagram of the algorithm used is presented in Fig. 2. The steps performed (that are basically the same ones used in the standard GS algorithm), are the following:

- 0 At the beginning, a field with an amplitude given by the square root of the expected irradiance and a constant phase is taken.
- 1. The field is propagated from the image-plane to the object-plane.
- 2. The amplitude information is discarded, leaving only the phase information (for the phase mask).
- 3. The amplitude and phase of the illumination field are added to the phase information to obtain the resulting object field.
- 4. The field is propagated from the object plane to the image plane.
- 5. The resulting reconstructed image (square of the field amplitude) is compared with the expected one. By using the correlation between both images as a criterion, a decision is taken to finish the process or continue iterating.
- 6. The phase from the reconstructed image is combined with the field amplitude obtained from the expected irradiance, and the process is repeated from step 1

## 4. Details of the system used in the numerical calculations

To get the numerical results the optical system shown in Fig. 3 is simulated by using the following components: a spatial filter (SF) (a point source), a focusing lens (FL) (200 mm focal-length doublet - Edmund Scientific, 45179), a 50% 50% beam-splitting cube (BS), a reflective phase mask (PM), a microscope objective (MO) (taken from US Patent 3893751 [10]), and an image plane (P) taken as a flat surface. The distances used were 487 mm from PM to the focusing point (FP) of FL, and 160 mm from FP to MO (to keep the correct working distance of MO).

The proposed system is highly aberrated as can be seen in Fig. 4(a). To test that it is also translational-variant, the first step was to calculate the phase mask needed to correct the aberrations and obtain a diffraction-limited system for an on-axis image [Fig. 4(b)]. Afterward, a tilt was applied to the phase mask to displace the image point along the X axis. As can be seen in Fig. 5, the spot radius changed when displaced, showing the system is translational-variant.

Then, to prove that the standard GS algorithm is not adequate for the optical system used, a phase mask using Fig. 6(a) as the target image was calculated. It can be seen in the reconstructed image [Fig. 6(b)], that the standard GS algorithm did not produce an adequate result in this case.

Finally, a test of the algorithm proposed in this work was conducted by using the same target image as in the previous tests [Fig. 6(a)]. In Fig. 7 some intermediate results of the iterations can be seen.

In Fig. 8, the normalized correlation between the reconstructed images and the target image as a function of the number of iterations of the algorithm is shown. It can be seen that most of the improvement in reconstruction happens, for this test, in the first 10 iterations and then the correlation stagnates.

## 5. Conclusions

In this work, it was shown that it is possible to use a GS-like algorithm to calculate a phase mask to be used in an aberrated and translational-variant optical system to produce an arbitrary image. The result was obtained by using a mixed algorithm, where the propagation of the optical field is calculated by raytracing point sources and using the Huygens-Fresnel principle. One assumption being made in this process is that the geometric phase can be regarded as the physical phase of the wave, and for this to be true, the Eikonal [11] equation must be valid. It can be proven that it is valid for plane and spherical waves [12] without regarding the wavelength, and approximately valid for Gaussian beams, and other slowly varying phase functions. For setups where the Eikonal equation is not valid on the image plane, the raytracing propagation can be performed up to a plane where the phase varies slowly enough for the Eikonal equation to be approximately valid, and then propagated back down to the observation plane using any of the standard methods to solve the diffraction integrals.

## References

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

**2. **J. S. Liu and M. R. Taghizadeh, “Iterative algorithm for the design of diffractive phase elements for laser beam shaping,” Opt. Lett. **27**, 1463–1465 (2002) [CrossRef] .

**3. **Q. Li, H. Gao, Y. Dong, Z. Shen, and Q. Wang, “Investigation of diffractive optical element for shaping a Gaussian beam into a ring-shaped pattern,” Opt. Laser Technol. **30**, 511–514 (1998) [CrossRef] .

**4. **C. Bay, N. Hubner, J. Freeman, and T. Wilkinson, “Maskless photolithography via holographic optical projection,” Opt. Lett. **35**, 2230–2232 (2010) [CrossRef] [PubMed] .

**5. **D. C. O’Shea, A. D. Kathman, and D. W. Prather, *Diffractive Optics: Design, Fabrication, and Test* (SPIE Publications, 2003) [CrossRef] .

**6. **O. Ripoll, K. Ville, and H. P. Herzig, “Review of iterative Fourier-transform algorithms for beam shaping applications,” Opt. Eng. **43**, 2549–2556 (2004) [CrossRef] .

**7. **G. Yang, B. Gu, J. Zhuang, and O. K. Ersoy, “Gerchberg-Saxton and Yang-Gu algorithms for phase retrieval in a nonunitary transform system: a comparison,” Appl. Opt. **33**, 209–218 (1994) [CrossRef] [PubMed] .

**8. **Z. Zalevsky and D. Mendlovic, “Gerchberg-Saxton algorithm applied in the fractional Fourier or the Fresnel domain,” Opt. Lett. **21**, 842–844 (1996) [CrossRef] [PubMed] .

**9. **J. W. Goodman, *Introduction to Fourier Optics* (Roberts & Co, 2005).

**10. **A. Shoemaker, 40X Microscope objective, US Patent 3893751 (1975).

**11. **M. Born, E. Wolf, and A. B. Bhatia, *Principles of Optics* (Pergamon Press, 1975).

**12. **N. Lindlein, “Simulation of micro-optical systems including microlens arrays,” Journal of Optics A: Pure and Applied Optics **4**, S1–S9 (2002) [CrossRef] .