## Abstract

We propose a novel design method for a circularly symmetric phase mask to extend the depth of focus. Using the free-form rational function as the solution space, we optimize the profile of the phase mask by analysis of the axial intensity distribution, which can be calculated efficiently by employing the fast Fourier transform algorithm. Numerical comparisons prove the resulting rational phase mask’s superiority to the existing quartic phase mask in intensity distribution and imaging performance.

© 2009 Optical Society of America

## 1. INTRODUCTION

Conventional imaging systems often suffer from defocus error, and various optical elements have been proposed to extend the depth of focus (DOF). The apodizer [1, 2], which is an amplitude mask, alleviates the defocus error at the expense of the light transmission and results in a decreased contrast. The odd symmetric phase masks can achieve a large focal depth after digital image restoration without energy loss [3, 4, 5]. However, this approach introduces a phase modulation changing continuously; the resulting artifacts make the image restoration difficult and time-consuming [6]. Also, this kind of odd symmetric phase mask is difficult and expensive to manufacture and test.

Circularly symmetric pure phase elements were also proposed in the literature, such as the logarithmic axicon (LA), the quartic phase mask (QPM), and the refractive and diffractive hybrid lens [7, 8, 9]. The last is fabricated through complicated laser writing on photoresist and subsequent reactive ion etching [9], while the former two can be more conveniently manufactured by use of traditional techniques such as diamond machining, and the cost can be reduced. Unlike the odd symmetric phase masks whose point-spread functions (PSF) are asymmetric, the PSFs of these circularly symmetric phase masks tend to be circularly symmetric as well. Additionally, the PSFs are also sharp and well localized for a large range of defocus without any shifts in the image sensor. Therefore, the post digital image restoration can be omitted and the circularly symmetric phase masks can satisfy some real-time all-optical applications to extend the DOF.

The LA was derived on the basis of geometrical optics: the geometrical law of energy conservation in [7] and Fermat’s principle in [10]. The QPM was achieved from the Wigner distribution function in [8] and the stationary phase method in [11]. They have similar performance in extending the DOF, and comparisons proved that the QPM is even superior to the LA [12]. However, there still exist two main problems with them. First, numerical simulations indicate that the Fresnel diffraction integral shows prominent rapid oscillations of the axial optical intensity, especially near the front and the end of the focal region, which are caused by diffraction and interference effects [13]. Second, near the front of the focal region the intensity distribution across a transverse plane is broad, monotonically decreasing away from the axis, and possesses no zeros. The two techniques also do not preserve uniformity of lateral resolution and energy flow along the focal range [13]. To solve these problems, the apodized annular-aperture approach has been proposed and widely used [14, 15], although it brings the drawbacks of a decrease of signal-to-noise ratio (SNR) and a possible decrease in image resolution. Some other possible approaches, such as employing partially coherent or polychromatic light, have been suggested and analyzed [16, 17], but strict illumination conditions limit their application [18].

In fact, there is the potential to improve these poor performances by optimizing the profile of the phase mask. In this paper we describe the optimization design method of a circularly symmetric phase mask called the rational phase mask (RPM) and compare it with the well-known QPM derived from analytical methods. We find this numerical approach provides more design freedom and allows obtaining improved performance with wide applications. Section 2 gives the basic theory. Using the free-form rational function as the solution space, we present the optimization design method of the RPM by analysis of the axial intensity distribution in Section 3. Section 4 is the numerical comparisons of intensity distribution and imaging performance between the optimized RPM and the QPM. We conclude the work in Section 5.

## 2. THEORY

We consider a circularly symmetric optical system with the generalized phase pupil function [19]

*ρ*is the radial coordinate at the exit pupil plane and $\varphi \left(\rho \right)$ represents the phase function. From the diffraction theory, the intensity PSF can be written as [19]

*r*is the radial coordinate at the image plane, ${J}_{0}(\mathbf{\cdot})$ denotes a Bessel function of the first kind and zero order, and the defocus parameter

*ψ*is defined as [19]

*L*is the aperture diameter;

*λ*is the wavelength; $f,{d}_{o},{d}_{i}$ are the focal length, the object distance, and the image sensor plane distance, respectively; and ${W}_{20}$ is the maximum phase shift. Then the axial intensity is

To extend the DOF, we usually expect a rectangular axial intensity distribution with a mean value as high as possible in the focal region. However, we observe that from the Parseval theorem [19], the area under the axial intensity distribution curve is

## 3. OPTIMIZATION DESIGN

We may take the work of pursuing the best phase mask as the approximation of the unknown optimal phase function $\varphi \left(\rho \right)$. Therefore the search space of the feasible solution should be as large as possible. Usually, a polynomial function is used to approximate an analytical function, but it often suffers from unsatisfactory convergence and a large number of undetermined coefficients, which means an increasing computational cost in the optimization procedure. The most advanced approximant would belong to a class of functions known as the rational function [20], which can obtain sufficient accuracy with a relatively small number of approximation coefficients and exhibits more flexible behavior than polynomials. Here we present the phase mask design method using the rational function, and the corresponding phase function for this RPM is described as [20]

Fig. 1 shows a typical axial intensity distribution curve resulting from Eq. (5). We assume the axial intensity oscillation is measured as the peak-to-valley value in the effective focal region and define a constant *C* as the ratio of the acceptable axial intensity oscillation to the maximum axial intensity value ${I}_{\mathrm{max}}$. The effective focal region $\Delta f$ can be defined as the region where there is no axial intensity value below $(1-C){I}_{\mathrm{max}}$. Then the merit function is to maximize $\Delta f$ under required constraints, which can be written as

*K*constrains the maximum axial intensity value to achieve a high enough SNR. For the convenience of comparison, the in-focus axial intensity value of a clear aperture (CA) is designated as unity. In this paper, $K=0.05$ $(C=0.1)$, 0.1 $(C=0.1)$, and 0.25 $(C=0.025)$ are considered and the corresponding masks are denoted by RPM0.05, RPM0.1, and RPM0.25, respectively. The simulated annealing algorithm is adopted to solve this constrained nonlinear optimization problem [21], and the resulting RPM phase profiles are plotted in Fig. 2 .

The phase function of the well-known QPM is given by [8, 11]

where the coefficient*α*is related to the DOF and a constant term is omitted. For comparison, we let $\alpha =10.95,5.45,2.36$ and achieve the comparable QPMs with the ${I}_{\mathit{\text{mean}}}$ equivalent to the optimized RPMs’. We denote them as QPM0.05, QPM0.1, and QPM0.25, respectively, and the corresponding phase profiles are also plotted in Fig. 2. As can be seen, all the slopes of these phase profiles are comparable, which implies these masks may have a similar manufacturability. In Section 4, we will compare them in detail.

## 4. NUMERICAL RESULTS

#### 4A. Intensity Distribution

From Eqs. (4, 2), the axial intensity distribution and the partial PSF near the optical axis for each of the three phase masks are calculated, respectively, as shown in Figs. 3, 4, 5 . Here, the radial coordinate *r* at the image plane is measured in units of the radius to the first zero of the CA Airy pattern for the in-focus diffraction-limited case, which is given by [19]

*α*. The axial intensity uniformity for the RPM0.05 and RPM0.1 is improved visibly, which results from the small value of

*C*. If we choose a larger

*C*(e.g., $C=0.5$), the optimized RPMs tend to converge to the corresponding QPMs. Particularly, for the QPM0.25 in Fig. 5a, the oscillations are relatively small, which reach about 13% of the maximum value. However, it is still possible to smooth the axial intensity distribution by optimizing the profile of the RPM. By setting $C=0.025$, we achieved the RPM0.25 with an extremely flat top in the focal region.

Simultaneously, the axial intensity uniformity is achieved at a cost of some decrease of the effective focal range compared with the QPM. This phenomenon was also observed when using the apodized annular-aperture approach [15], where ${\Delta}_{z}$ was defined as the transition region between the top and the bottom of the axial intensity distribution curve, and a larger ${\Delta}_{z}$ resulted in a better smoothness. In Fig. 5a we also plot the axial intensity distribution of the CA. The extending of the DOF is obvious and Eq. (6) is verified, that is, we cannot achieve a large DOF and a high SNR at the same time.

We also evaluate the off-axis quality of the light. From the partial PSF plotted in Figs. 3, 4, the QPMs’ PSF side lobes or local maxima show some discontinuity along the focal region. Usually a larger central peak corresponds to a smaller second maximum, and a larger third maximum, and so on. The similarity between the PSFs at different image planes is destroyed. For the partial PSF of the RPMs, the local maxima are more continuous along the focal region and the lateral intensity uniformity is preserved better. But for the RPM0.25 in Fig. 5b, this improvement is not so evident. The description above may hold true except for the front of the focal region where the central spot size increases and the local maxima disappear. In the paraxial region of the partial PSF, we calculate the light energy,

that flows through a circular hole with the constant radius ${r}_{m}$ for the RPM0.1 and QPM0.1. If ${r}_{m}$ is reduced to near zero, we can obtain an approximation of the axial intensity. Here, we let ${r}_{m}=0.64{r}_{0}$, which is the minimum radius from the axis to the first zero intensity along the focal region. The resulting energy flow curves are plotted as the solid ones in Fig. 6 . For the QPM0.1, the nonuniform character of the energy flow along the focal region consists of oscillations and a decline. For the RPM0.1 the main trend is only the decline, and the value between the maximum and the minimum is smaller. The dot curves in Fig. 6 correspond to the results calculated from letting ${r}_{m}=0.32{r}_{0}$, and the uniformity is improved significantly for both RPM0.1 and QPM0.1. This uniformity of the energy flow may be quite useful in applications such as energy concentration, laser machining, and so on.#### 4B. Imaging Performance

We have considered the axial and paraxial irradiance, but this common analysis is incomplete for imaging applications that require full information about the PSF or the equivalent optical transfer function (OTF). Next, taking RPM0.1 and QPM0.1, for instance, we concern ourselves with the radius to the first zero of the PSF, which we label as the cutoff lateral position ${r}_{c}$, and we may expect that ${r}_{c}$ would be approximately constant along the focal region. However, Fig. 7a illustrates that near the front of the focal region $(\psi <-32)$ the PSF profiles of QPM0.1 are apt to decrease monotonically and possess no zeros, as pointed out in [13], while for the RPM0.1, its PSF profiles always possess zeros. Elsewhere, the cutoff lateral positions of the RPM0.1 and QPM0.1 tend to converge to the same value ${r}_{c}=0.64{r}_{0}$ with tiny fluctuations. It is interesting to point out that ${r}_{c}<{r}_{0}$ does not mean the RPM0.1 or QPM0.1 has a higher resolution than the in-focus CA, because the local maxima contain energy no less than the central spot as shown in Fig. 7b. For instance, at $\psi =35$ the central peak contains only about 1% of the total energy for both the RPM0.1 and QPM0.1; for comparison, at $\psi =-35$ the central peak contains about 16% of the total energy for both the RPM0.1 and QPM0.1. This can also be seen from the modulation transfer functions (MTF) plotted in Fig. 8 , where the cutoff frequencies are very low for $\psi >0$. For $-10<\psi <0$ and $-60<\psi <-30$, the cutoff frequencies of the RPM0.1 are higher. And for $-30<\psi <-10$, the cutoff frequencies are almost the same for the two masks. It is notable that all the cutoff frequencies are reduced compared with the those of the in-focus CAs, as shown in Fig. 10a. For applications using pixelated detectors, these masks may be usefully employed as antialiasing filters [22].

We also compare the simulated imaging performance of the two masks. Figure 9 shows the raw images of a simulated spoke target produced by applying the corresponding OTFs, and there is no other digital image processing. At $\psi =-55$, $-45$, and $-35$, the imaging performance of the RPM0.1 is better than that of the QPM0.1, and more details can be seen. At $\psi =-25$ and $-15$, the results are comparable, while at $\psi =-5$ the high-frequency components in the center disappear first for the QPM0.1, and at $\psi =5$ the image degradation becomes serious for both masks, and most details are lost. For comparison, the simulated images of the CA are given in Fig. 10b and the high-frequency cutoff and contrast reversal become evident at $\psi =\pm 5$ and $\pm 7$. The RPM0.1 has extended the DOF to a larger range than the QPM0.1, but they both have a contrast decrease resulting from their low MTFs in the focal region. The simulation above is under the noise-free condition, and in future research, the effects of noise should be taken into consideration.

## 5. CONCLUSION

In conclusion, we have proposed a novel design method of the circularly symmetric phase mask to extend the DOF by analysis of the axial intensity distribution, which can be calculated efficiently by employing the fast Fourier transform algorithm. The simple constraints on the axial intensity oscillations result in improved uniformity of the axial intensity, energy flow, lateral resolution, and also improved imaging performance along the focal region. With the free-form rational function as the solution space, the numerical optimization method may present more design flexibility for various applications by more complex statistical analysis of the axial intensity distribution.

**1. **J. Ojeda-Castenada, P. Andres, and A. Diaz, “Annular apodizers for low sensitivity to defocus and to spherical aberration,” Opt. Lett. **11**, 487–489 (1986). [CrossRef]

**2. **J. Ojeda-Castenada, E. Tepichin, and A. Diaz, “Arbitrary high focal depth with a quasi-optimum real and positive transmittance apodizer,” Appl. Opt. **28**, 2666–2670 (1989). [CrossRef]

**3. **E. R. Dowski and W. T. Cathey, “Extended depth of field through wave-front coding,” Appl. Opt. **34**, 1859–1866 (1995). [CrossRef] [PubMed]

**4. **Q. Yang, L. Liu, and J. Sun, “Optimized phase pupil masks for extended depth of field,” Opt. Commun. **272**, 56–66 (2007). [CrossRef]

**5. **S. Bagheri, P. Silveira, and D. Farias, “Analytical optical solution of the extension of the depth of field using cubic-phase wavefront coding. Part II. Design and optimization of the cubic phase,” J. Opt. Soc. Am. A **25**, 1064–1074 (2008). [CrossRef]

**6. **S. Förster, H. Gross, F. Höller, and L. Höring, “Extended depth of focus as a process of pupil manipulation,” Proc. SPIE **5962**, 596207 (2005). [CrossRef]

**7. **J. Sochacki, S. Bara, Z. Jaroszewicz, and A. Kołodziejczyk, “Phase retardation of uniform-intensity axilens,” Opt. Lett. **17**, 7–9 (1992). [CrossRef] [PubMed]

**8. **D. Zalvidea and E. E. Sicre, “Phase pupil functions for focal depth enhancement derived from a Wigner distribution function,” Appl. Opt. **37**, 3623–3627 (1998). [CrossRef]

**9. **Z. Liu, A. Flores, M. R. Wang, and J. J. Yang, “Diffractive infrared lens with extended depth of focus,” Opt. Eng. **46**, 018002 (2007). [CrossRef]

**10. **W. Chi and N. George, “Electric imaging using a logarithmic asphere,” Opt. Lett. **26**, 875–877 (2001). [CrossRef]

**11. **S. Mezouari and A. R. Harvey, “Phase pupil functions for control defocus and spherical aberrations,” Opt. Lett. **28**, 771–773 (2003). [CrossRef] [PubMed]

**12. **D. Zalvidea, C. Colautti, and E. E. Sicre, “Quality parameters analysis of optical imaging systems with enhanced focal depth using the Wigner distribution function,” J. Opt. Soc. Am. A **17**, 867–873 (2000). [CrossRef]

**13. **L. R. Staronski, J. Sochacki, Z. Jaroszewicz, and A. Kołodziejczyk, “Lateral distribution and flow of energy in uniform-intensity axicons,” J. Opt. Soc. Am. A **9**, 2091–2094 (1992). [CrossRef]

**14. **J. Sochacki, Z. Jaroszewicz, L. R. Staronski, and A. Kołodziejczyk, “Annular-aperture logarithmic axicons,” J. Opt. Soc. Am. A **10**, 1765–1768 (1993). [CrossRef]

**15. **Z. Jaroszewicz, J. Sochacki, A. Kołodziejczyk, and L. R. Staronski, “Apodized annular-aperture logarithmic axicon: smoothness and uniformity of intensity distributions,” Opt. Lett. **18**, 1893–1895 (1993). [CrossRef] [PubMed]

**16. **S. Y. Popov and A. T. Friberg “Design of diffractive axicons for partially coherent light,” Opt. Lett. **23**, 1639–1641 (1998). [CrossRef]

**17. **Z. Jaroszewicz, J. F. Roman Dopazo, and C. Gomez Reino, “Uniformization of the axial intensity of diffraction axicons by polychromatic illumination,” Appl. Opt. **35**, 1025–1031 (1996). [CrossRef] [PubMed]

**18. **S. Y. Popov and A. T. Friberg, “Apodization of generalized axicons to produce uniform axial line images,” Pure Appl. Opt. **7**, 537–548 (1998). [CrossRef]

**19. **J. W. Goodman, *Introduction to Fourier Optics*, 3rd ed. (Roberts, 2005).

**20. **J. H. Mathews and K. D. Fink, *Numerical Methods Using MATLAB*, 4th ed. (Prentice Hall, 2004).

**21. **S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science **220**, 671–680 (1983). [CrossRef] [PubMed]

**22. **S. Mezouari, G. Muyo, and A. R. Harvey, “Circularly symmetric phase filters for control of primary third-order aberrations: coma and astigmatism,” J. Opt. Soc. Am. A **23**, 1058–1062 (2006). [CrossRef]