1. INTRODUCTION

Wavefront coding is an optical–digital imaging technique first proposed by Dowski and Cathey [1]. A standard fixed-focus imaging system combined with coding optics, e.g., a phase filter in the pupil plane, produces an encoded image (on a photosensor), which is subsequently decoded by digital processing. The result of such a coding–decoding procedure is a final image with extended depth of focus [1, 2, 3, 4, 5]. Similar approaches based on modifying the phase of the pupil function to make the imaging system invariant to specific wavefront aberrations, including defocus, are reported in [6, 7, 8, 9]. Clearly, the defocus aberration is of major interest because of its high potential for technical and consumer imaging applications.

Analytical expressions for the point spread function or, alternatively, the optical transfer function (OTF) of an imaging system with an encoding phase mask enable characterization of the system behavior in the physical domain (in terms of the point spread function) or, alternatively, in the spatial frequency domain (in terms of the OTF). Moreover, such an expression permits evaluation of the optimal parameters of the phase mask by, for example, maximizing the modulation transfer function (MTF) in the spatial frequency range of interest. Most theoretical studies of wavefront coding assume a rectangular aperture (in the two-dimensional case) and mathematical separability of the pupil function [5, 10, 11, 12, 13]. This allows the OTF, specified by the double integral over the pupil region, to be reduced to a product of two single integrals, which, for example, are approximately evaluated by using the stationary phase method as suggested in [1]. However, additional assumptions are needed to derive analytical expressions for optical systems with circular or other shaped apertures and/or for systems with a nonseparable pupil function, since in this case the two-dimensional integrals cannot be evaluated as a product of two one-variable integrals [14].

In this paper, we propose a two-dimensional generalization of the methods developed for approximating the OTF of a one-dimensional fixed-focus optical system with an encoding phase mask. The analysis is based on the multidimensional stationary phase method, as described, for example, in [15, 16, 17], and the analysis does not require separability of the pupil function. An approximate expression for the OTF of the imaging system invariant to second-order aberrations resulting from a third-order phase mask is obtained explicitly. To the best of our knowledge, this is the first study of the two-dimensional OTF asymptotics that does not require additional simplifications and assumptions about the pupil function and the pupil shape.

2. GENERAL ANALYSIS

Consider an optical incoherent system with a phase mask in the pupil plane specified by an unknown function $\theta (x,y)$ in the Hopkins’ canonical pupil coordinates x and y [18]. Following [1], we look for the function $\theta (x,y)$ that makes the system OTF invariant with respect to quadratic phase distortion represented by

In accordance with [18, 19], the two-dimensional OTF of an incoherent imaging system, as a function of the magnitude $w_0$ of the quadratic phase distortion, can be expressed as a normalized autocorrelation function of the complex pupil function

The aim of our analysis is to find such a phase mask function $\theta (x,y)$ that, first, makes the system OTF $H({\omega }_x,{\omega }_y,w_0)$ weakly dependent on $w_0$ over a wide range of $w_0$ and, second, ensures that $H({\omega }_x,{\omega }_y,w_0)$ can be inverted without excessive noise amplification at all spatial frequencies ${\omega }_x$, ${\omega }_y$.

To simplify the analysis of the OTF of the imaging system, we introduce the following notation:

Assuming that the function $F(x,y,{\omega }_x,{\omega }_y)$ changes rapidly with x and y, the integral Eq. (9) can be evaluated approximately by using the two-dimensional stationary phase method [15, 16]. Leaving the dominant term in the approximation, Eq. (9) yields

The quadratic form generated by the symmetric real matrix M specified by Eq. (12) can be converted to a diagonal form by a nondegenerate transformation. In conformity with Sylvester’s law of inertia [20], the difference between the number of positive and negative coefficients in the diagonal representation of the matrix, i.e., the matrix signature, remains constant and does not depend on the transformation. In particular, the quadratic form generated by Eq. (12) can be diagonalized as follows: ${\lambda }_1{\tilde{x}}^2+{\lambda }_2{\tilde{y}}^2$, where ${\lambda }_1$ and ${\lambda }_2$ are the real eigenvalues of the symmetric matrix M and $\tilde{x}$, $\tilde{y}$ are the pupil coordinates in a new orthogonal basis [20].

If the phase function $\Phi =F(x,y,{\omega }_x,{\omega }_y)+2\pi (x{\upsilon }_x+y{\upsilon }_y)$ has several stationary points, then the additional terms similar to those on the right-hand side of Eq. (10) should be added to the OTF approximation.

It is worth noting that the transformation $L:\{(x,y),F\}\to \{({\upsilon }_x,{\upsilon }_y),\tilde{\Phi }\}$, where $F=F(x,y,{\omega }_x,{\omega }_y)$, $\tilde{\Phi }=F(x_c,y_c,{\omega }_x,{\omega }_y)+2\pi (x_c{\upsilon }_x+y_c{\upsilon }_y)$, and $x_c$, $y_c$ are defined in terms of ${\upsilon }_x$, ${\upsilon }_y$ according to Eqs. (11), is per se the Legendre transformation [16].

From Eqs. (10, 11, 12, 13) it is clear that the OTF modulus $\mid \tilde{H}({\omega }_x,{\omega }_y,{\upsilon }_x,{\upsilon }_y)\mid$, i.e., the MTF of the imaging system, is specified only by the Hessian $\Delta (x_c,y_c)$. We choose the phase mask function $\theta (x,y)$ in such a way that $\Delta (x_c,y_c)$ does not depend on ${\upsilon }_x$, ${\upsilon }_y$ and, hence, on the magnitude $w_0$ of the quadratic phase distortion. Thus, $\Delta (x_c,y_c)$ satisfies the conditions

3. CUBIC PHASE MASK

Consider now the phase mask function $\theta (x,y)$ chosen in the form of a third-order polynomial

For lower-order polynomials, for example $\theta (x,y)=\alpha xy+\beta (x^2+y^2)$, from Eqs. (11) it follows that there are no stationary points for the function $F(x,y,{\omega }_x,{\omega }_y)$. For higher-order polynomials, the stationary points are described by a system of two nonlinear equations (11), and the dependences of $x_c({\upsilon }_x,{\upsilon }_y)$ and $y_c({\upsilon }_x,{\upsilon }_y)$ on the parameters ${\upsilon }_x$, ${\upsilon }_y$ become nonlinear; i.e., the conditions specified by Eqs. (20) are no longer fulfilled (see Appendix A for details). Thus, it can be concluded that the OTF approximated by Eq. (10) does not depend on the magnitude $w_0$ of the quadratic phase distortion only when the phase mask $\theta (x,y)$ is represented by a third-order polynomial in the pupil coordinates x and y.

For the third-order phase mask specified by Eq. (21), Eq. (11) results in

In the case of $\alpha \to \infty$, from Eqs. (24, 25) it follows that the stationary point $(x_c,y_c)$ tends to the origin ($x=0$, $y=0$), and, thus, the stationary point lies within the pupil area Ω. From Eq. (22), the Hessian matrix M, defined by Eq. (12), at $(x_c,y_c)$ becomes

Finally, substituting Eqs. (8, 22, 24, 25, 26) into Eq. (10), the approximation of the two-dimensional OTF of the imaging system invariant to second-order phase distortions (at large values of α and ${\omega }_x^2+{\omega }_y^2\ne 0$) takes the form

The precision of the OTF asymptotic representation, Eq. (30), at large α is determined by the residual term ${\mid \Delta \mid }^{-1∕2}O\left({\alpha }^{-1}\right)$ or $O\left({\alpha }^{-2}\right)$. However, a higher precision can be achieved by evaluating additional terms of the order $O\left({\alpha }^{-N}\right)$, where $N⩾2$ is an integer [15, 16]. Alternatively, the contribution from the exponentially decaying terms in oscillatory integrals can be used to improve the accuracy of the OTF approximation, Eq. (10) [21]. This approach provides reasonable accuracy even for moderate values of the asymptotic parameter, in our case α.

Using the general formula, Eq. (31), which approximates the OTF of the optical system with the generalized cubic mask, we now consider two particular examples: (i) a phase mask with mixed x and y terms, ${\gamma }_1=0$, ${\gamma }_2=0$, and (ii) a cubic phase mask with no mixed terms, ${\beta }_1={\beta }_2=0$, ${\gamma }_1={\gamma }_2=1$.

In the first case, (i), Eq. (21) simplifies to $\theta (x,y)=\alpha ({\beta }_1{x}^{2}y+{\beta }_2{y}^{2}x)$, and the phase function contains only mixed x and y terms. In agreement with Eq. (27), the denominators of Eq. (24, 25, 26) are nonzero at all values of ${\beta }_1$, ${\beta }_2$ and at ${\omega }_x^2+{\omega }_y^2\ne 0$. One may conclude from Eq. (29) that the eigenvalues ${\lambda }_1$, ${\lambda }_2$ of the Hessian matrix M have different signs at all ${\beta }_1$, ${\beta }_2$ and, thus, $\delta (x_c,y_c)=0$. Analogously to Eq. (31), the OTF becomes

In the second case, (ii), the phase function is $\theta (x,y)=\alpha (x^3+y^3)$, which corresponds to the cubic phase mask originally suggested in [1]. From Eq. (29) ${\lambda }_1=3{\omega }_x$, ${\lambda }_2=3{\omega }_y$, and the signature of M evaluates to $\delta =\mathrm{sgn}\left({\omega }_x\right)+\mathrm{sgn}\left({\omega }_y\right)$. Thus, Eq. (30) at $\alpha \to \infty$ and ${\omega }_x\ne 0$, ${\omega }_y\ne 0$ gives

4. DISCUSSION AND COMPARISON WITH DIRECT NUMERICAL SIMULATIONS

Asymptotic equation (10) presents the two-dimensional OTF $H({\omega }_x,{\omega }_y,w_0)$ of the incoherent fixed-focus imaging system subject to the quadratic phase distortion $W(x,y)$, given by Eq. (1), and having a phase-coding mask in its pupil plane. With the phase $\theta (x,y)$ of the mask expressed by a third-order polynomial, Eq. (21), the OTF asymptotics at $\alpha \to \infty$ result in Eq. (30), which is invariant with respect to the magnitude $w_0$ of the phase distortion $W(x,y)$. The explicit analytical expression, Eq. (31), facilitates the evaluation of response properties of the incoherent optical system with the cubic mask. So, the MTF of a particular optical system can be easily optimized (according to requirements) by a proper choice of the parameters ${\beta }_1$, ${\beta }_2$, ${\gamma }_1$, and ${\gamma }_2$ of the mask.

The derivation of Eqs. (10, 30) implies that the phase function $\theta (x,y)$ varies rapidly with x and y, for example, because of large α in Eq. (21), and no other assumptions regarding the pupil shape and separability of the pupil function are used. As is seen from Eq. (30), at ${\beta }_1\ne 0$, ${\beta }_2\ne 0$ and with the parameters of the mask satisfying Eq. (27), the OTF approximation is a smooth function with no singularities at spatial frequencies ${\omega }_x^2+{\omega }_y^2\ne 0$. In contrast to Eq. (30), Eq. (33) together with Eqs. (34, 35) exhibit singularities at ${\omega }_x=0$ and ${\omega }_y=0$. In these regions, the OTF requires additional normalization to unity.

To validate the stationary-phase approximation of the system OTF and assess the accuracy of the method, a comparison has been made between the MTFs, i.e., $\mid H({\omega }_x,{\omega }_y,w_0)\mid$, evaluated with Eq. (31) and purely numerical calculations based on Fourier transformation of the pupil function $P(x,y)$ with the cubic phase mask $\theta (x,y)$, defined by Eqs. (3, 21), respectively. The numerical calculations have been performed by using the OTF expressed via $P(x,y)$, for example [19],

Figure 1 shows the MTF cross sections at ${\omega }_y=0$ calculated according to Eqs. (31, 36) for phase mask parameters $\alpha =10,20,50$ and varying degree of defocus $w_0=0,5,10,20$. Simulations have been carried out for three phase masks: plots 1–3, $\theta (x,y)=\alpha (x^3+y^3)$; plots 4–6, $\theta (x,y)=\alpha (x^{2}y+y^{2}x)$; plots 7–9, $\theta (x,y)=\alpha (x^{2}y+y^{2}x+x^3+y^3)$. In conformity with the stationary phase method, which requires $\alpha \to \infty$, plots at $\alpha =50$ show better agreement with the asymptotic curves than plots corresponding to $\alpha =10,20$. It can be also seen that at frequencies $1⩽\mid {\omega }_x\mid ⩽2$ the phase functions with mixed terms $x^{2}y$ and $x{y}^2$ demonstrate better concordance with the approximations. One may conclude that a cubic phase filter with properly chosen $x^{2}y$ and $x{y}^2$ terms may eventually require a lower signal-to-noise premium in comparison with the phase mask function $\theta (x,y)=\alpha (x^3+y^3)$. We should note that an enhanced signal-to-noise ratio and superior imaging quality with the mixed-term phase mask were reported earlier in [7].

From Eq. (31), a restoration scheme corresponding, for example, to the least-mean-square-error filter can obviously be found [19]. The spatial spectrum of an object is then calculated from the intermediate image encoded with the phase mask $\theta (x,y)$.

Often the phase mask function $\theta (x,y)$ is expressed in polynomial form in the canonical pupil coordinates x and y [7]. From the above theoretical analysis (see Appendix A for details) it follows that the OTF approximation does not depend on the magnitude $w_0$ of the quadratic phase distortion, defined by Eq. (1), only if the phase mask $\theta (x,y)$ is a third-order polynomial in x and y. For lower- and higher-order polynomials either the system of equations (11) or conditions (20) are not fulfilled.

Note that approximate solutions in the form of logarithmic [4] and odd-symmetric quadratic [11] phase functions have recently been reported for one-dimensional (mathematically separable) phase masks. Analysis of such nonlinear solutions is beyond the scope of this study.

5. CONCLUSION

We have derived an approximate analytical expression for the two-dimensional OTF of a fixed-focus imaging system subject to quadratic phase distortions and having a phase mask in its pupil plane. For the phase mask function in polynomial form, the system response becomes invariant to quadratic phase distortions only if the phase is a third-order polynomial in the pupil coordinates. The OTF approximation is found to be a smooth function with no singularities at all spatial frequencies except zero, but only when the pupil phase function contains mixed terms in the pupil coordinates.

The approximate analytical expression can be used for computer simulation and design of phase filters for wavefront coding–decoding imaging systems, design of ophthalmic lenses providing extended depth of focus, and likely other technical and medical applications.

APPENDIX A: ENCODING FUNCTION FOR QUADRATIC PHASE DISTORTION

In this appendix we prove that for the quadratic phase distortion $W(x,y)$ specified by Eq. (1) the asymptotic representation of the OTF $H({\omega }_x,{\omega }_y,w_0)$ according to Eq. (10) becomes invariant with respect to the amplitude $w_0$ of $W(x,y)$ only if the phase mask function $\theta (x,y)$ is a third-order polynomial in x and y.

Invariance of $H({\omega }_x,{\omega }_y,w_0)$ on $w_0$ requires Eqs. (14) to be fulfilled. By recalling the Hessian definition, Eq. (13), Eqs. (14) can be rewritten in terms of $F(x,y,{\omega }_x,{\omega }_y)$, where $x=x_c({\upsilon }_x,{\upsilon }_y)$ and $y=y_c({\upsilon }_x,{\upsilon }_y)$ obey Eqs. (11),

(A1)

$$\frac{\partial }{\partial {\upsilon }_x}[\frac{{\partial }^{2}F}{\partial x^2}\frac{{\partial }^{2}F}{\partial y^2}-\frac{{\partial }^{2}F}{\partial x\partial y}\frac{{\partial }^{2}F}{\partial y\partial x}]=0,$$

$$\frac{\partial }{\partial {\upsilon }_y}[\frac{{\partial }^{2}F}{\partial x^2}\frac{{\partial }^{2}F}{\partial y^2}-\frac{{\partial }^{2}F}{\partial x\partial y}\frac{{\partial }^{2}F}{\partial y\partial x}]=0,$$

or, after differentiation,

(A2)

$$\frac{\partial x_c}{\partial {\upsilon }_x}[\frac{{\partial }^{3}F}{\partial x^3}\frac{{\partial }^{2}F}{\partial y^2}+\frac{{\partial }^{2}F}{\partial x^2}\frac{{\partial }^{3}F}{\partial y^2\partial x}-2\frac{{\partial }^{2}F}{\partial y\partial x}\frac{{\partial }^{3}F}{\partial x^2\partial y}]+\frac{\partial y_c}{\partial {\upsilon }_x}[\frac{{\partial }^{3}F}{\partial x^2\partial y}\frac{{\partial }^{2}F}{\partial y^2}+\frac{{\partial }^{2}F}{\partial x^2}\frac{{\partial }^{3}F}{\partial y^3}-2\frac{{\partial }^{2}F}{\partial x\partial y}\frac{{\partial }^{3}F}{\partial x\partial y^2}]=0,$$

$$\frac{\partial x_c}{\partial {\upsilon }_y}[\frac{{\partial }^{3}F}{\partial x^3}\frac{{\partial }^{2}F}{\partial y^2}+\frac{{\partial }^{2}F}{\partial x^2}\frac{{\partial }^{3}F}{\partial y^2\partial x}-2\frac{{\partial }^{2}F}{\partial y\partial x}\frac{{\partial }^{3}F}{\partial x^2\partial y}]+\frac{\partial y_c}{\partial {\upsilon }_y}[\frac{{\partial }^{3}F}{\partial x^2\partial y}\frac{{\partial }^{2}F}{\partial y^2}+\frac{{\partial }^{2}F}{\partial x^2}\frac{{\partial }^{3}F}{\partial y^3}-2\frac{{\partial }^{2}F}{\partial x\partial y}\frac{{\partial }^{3}F}{\partial x\partial y^2}]=0.$$

Equations (A2) are fulfilled if $x_c({\upsilon }_x,{\upsilon }_y)=\text{constant}$ and $x_y({\upsilon }_x,{\upsilon }_y)=\text{constant}$, or, alternatively, if the third- or second-order partial derivatives of F are zeros. In the first case, Eqs. (17) become inconsistent. In the second case, taking into account that $\Delta (x_c,y_c)\ne 0$ and, thus, second-order derivatives of F cannot be zeros, we obtain

(A3)$$F(x,y,{\omega }_x,{\omega }_y)=\sum _{n+m⩽2}{a}_{nm}({\omega }_x,{\omega }_y)x^n{y}^m,$$

where $a_{nm}=a_{nm}({\omega }_x,{\omega }_y)$ are the coefficients. On the other hand, by substituting the Taylor expansions of $\theta (x+{\omega }_x∕2,y+{\omega }_y∕2)$ and $\theta (x-{\omega }_x∕2,y-{\omega }_y∕2)$ at $x,y$ into Eq. (7), it can be found that

(A4)$$F(x,y,{\omega }_x,{\omega }_y)={\omega }_x\frac{\partial \theta (x,y)}{\partial x}+{\omega }_y\frac{\partial \theta (x,y)}{\partial y}+\frac{2}{3\text{!}}{(\frac{{\omega }_x}{2}\frac{\partial }{\partial x}+\frac{{\omega }_y}{2}\frac{\partial }{\partial y})}^3\theta (x,y)+\cdots .$$

By equating Eqs. (A3, A4) and performing integration over x and y, we obtain the phase mask function $\theta (x,y)$ in the form of a third-order polynomial (in x and y):

(A5)$$\theta (x,y)=\sum _{n+m⩽3}{b}_{nm}{x}^n{y}^m,$$

$b_{nm}$ being constants. The coefficients $a_{nm}$ in Eq. (A3) can be expressed in terms of $b_{nm}$,

(A6)

$$a_{00}=b_{12}{\omega }_x{\omega }_y^2+b_{21}{\omega }_x^2{\omega }_y∕4+b_{01}{\omega }_y+b_{10}{\omega }_x,$$

$$a_{01}=b_{11}{\omega }_x,$$

$$a_{10}=b_{11}{\omega }_y,$$

$$a_{11}=2b_{12}{\omega }_y+2b_{21}{\omega }_x,$$

$$a_{20}=2b_{21}{\omega }_y,$$

$$a_{02}=2b_{12}{\omega }_x.$$

Finally, the linear expressions $x=x_c({\upsilon }_x,{\upsilon }_y)$ and $y=y_c({\upsilon }_x,{\upsilon }_y)$ corresponding to Eq. (20) can be obtained by substituting Eq. (A3) into Eq. (11).

 

Fig. 1 MTF of the defocused optical system with the generalized cubic mask. Solid curves, MTF curves obtained with the approximate expression, Eq. (31); scatter plots represent the MTFs calculated numerically using the pupil function, as given by Eq. (36), for varying degrees of defocus $w_0$.

Download Full Size | PDF

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

2. J. van der Gracht, E. R. Dowski Jr., M. G. Taylor, and D. M. Deaver, “Broadband behavior of an optical–digital focus-invariant system,” Opt. Lett. 21, 919–921 (1996). [CrossRef]   [PubMed]  

3. S. Bradburn, W. T. Cathey, and E. R. Dowski, “Realization of focus invariance in optical–digital systems with wave-front coding,” Appl. Opt. 36, 9157–9166 (1997). [CrossRef]  

4. W. T. Cathey and E. R. Dowski, “New paradigm for imaging systems,” Appl. Opt. 41, 6080–6092 (2002). [CrossRef]   [PubMed]  

5. S. S. Sherif, W. T. Cathey, and E. R. Dowski, “Phase plate to extend the depth of field of incoherent hybrid imaging systems,” Appl. Opt. 43, 2709–2721 (2004). [CrossRef]   [PubMed]  

6. S. Mezouari and A. R. Harvey, “Phase pupil functions for reduction of defocus and spheric aberrations,” Opt. Lett. 28, 771–773 (2003). [CrossRef]   [PubMed]  

7. S. Prasad, T. C. Torgersen, V. P. Pauca, R. J. Plemmons, and J. van der Gracht, “High-resolution imaging using integrated optical systems,” Int. J. Imaging Syst. Technol. 14, 67–74 (2004). [CrossRef]  

8. A. Castro and J. Ojeda-Castaneda, “Asymmetric phase masks for extended depth of field,” Appl. Opt. 43, 3474–3479 (2004). [CrossRef]   [PubMed]  

9. 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]  

10. H. Lei, H. Feng, X. Tao, and Zh. Xu, “Imaging characteristics of a wavefront coding system with off-axis aberrations,” Appl. Opt. 45, 7255–7263 (2006). [CrossRef]   [PubMed]  

11. M. Somayaji and M. P. Christensen, “Frequency analysis of the wave-front coding odd-symmetric quadratic phase mask,” Appl. Opt. 46, 216–226 (2007). [CrossRef]   [PubMed]  

12. M. Somayaji and M. P. Christensen, “Enhancing form factor and light collection of multiplex imaging systems by using a cubic phase mask,” Appl. Opt. 45, 2911–2923 (2006). [CrossRef]   [PubMed]  

13. G. Muyo and A. R. Harvey, “Decomposition of the optical transfer function: wavefront coding imaging systems,” Opt. Lett. 30, 2715–2717 (2005). [CrossRef]   [PubMed]  

14. S. Bagheri, P. E. X. Silveria, and D. P. Farias, “Analytical optimal solution of the extension of the depth of filed using cubic-phase wavefront coding. Part I. Reduced-complexity approximate representation of the modulation transfer function,” J. Opt. Soc. Am. A 25, 1051–1063 (2008). [CrossRef]  

15. M. Fedoruk, Analysis I. Asymptotic Methods in Analysis, Vol. 13 of Encyclopaedia of Mathematical Science (Springer, 1989).

16. M. V. Fedoruk, Saddle Point Method (Nauka, 1977) (in Russian).

17. R. Wong, Asymptotic Approximations of Integrals (Academic, 1989).

18. H. H. Hopkins, “The frequency response of a defocused optical system,” Proc. R. Soc. London, Ser. A 231, 91–103 (1955).

19. J. W. Goodman, Introduction to Fourier Optics (Roberts & Company, 2005).

20. F. R. Gantmacher, The Theory of Matrices (Chelsea, 1959), Vol. 1.

21. D. Kaminski, “Exponentially improved stationary phase approximations for double integrals,” Methods Appl. Anal. 1, 44–56 (1994).

22. M. Frigo and S. G. Johnson, “FFTW,” http://www.fftw.org.