Abstract
An approximate analytical expression is derived for the two-dimensional incoherent optical transfer function (OTF) of an imaging system invariant to second-order aberrations. The system broadband behavior resulting from a third-order phase mask in its pupil plane is analyzed by using the two-dimensional stationary phase method. This approach does not require mathematical separability of the pupil function and can be applied to any pupil shape. The OTF is found to be a well-defined and smooth function at all nonzero spatial frequencies when the phase mask function includes third-order mixed terms in the pupil coordinates.
© 2009 Optical Society of America
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 in the Hopkins’ canonical pupil coordinates x and y [18]. Following [1], we look for the function that makes the system OTF invariant with respect to quadratic phase distortion represented by
in which the coefficient denotes the magnitude of the phase distortion and a, b, c are the constants satisfying , , . In particular, represents defocus when and , and represents astigmatisms when , , or when , .In accordance with [18, 19], the two-dimensional OTF of an incoherent imaging system, as a function of the magnitude of the quadratic phase distortion, can be expressed as a normalized autocorrelation function of the complex pupil function
where Ω is the pupil area, and are the reduced spatial frequencies [18], and the asterisk denotes complex conjugation. In Eq. (2) the generalized pupil function is given byHere is the imaginary unit, andBy substituting Eq. (3) into Eq. (2), one may obtainwherespecifies the area of overlap in the integral of Eq. (5), andis the phase contribution.The aim of our analysis is to find such a phase mask function that, first, makes the system OTF weakly dependent on over a wide range of and, second, ensures that can be inverted without excessive noise amplification at all spatial frequencies , .
To simplify the analysis of the OTF of the imaging system, we introduce the following notation:
and we rewrite the OTF defined by Eq. (5) in terms of , ,which is a Fourier transformation of the function .Assuming that the function 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
where (, ) is the extremum point, or alternatively, the stationary point, of the phase function , i.e., is derived fromThe Hessian in Eq. (10), i.e., the determinant of the Hessian matrixevaluated at (, ) is expressed asand is the signature, i.e., the difference between the number of positive and negative eigenvalues, of the Hessian matrix M at . As is seen from Eq. (10), the Hessian completely defines the behavior of , which is the system MTF. So, the mathematical analysis of will allow us to find a proper phase mask.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: , where and are the real eigenvalues of the symmetric matrix M and , are the pupil coordinates in a new orthogonal basis [20].
If the phase function 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 , where , , and , are defined in terms of , according to Eqs. (11), is per se the Legendre transformation [16].
From Eqs. (10, 11, 12, 13) it is clear that the OTF modulus , i.e., the MTF of the imaging system, is specified only by the Hessian . We choose the phase mask function in such a way that does not depend on , and, hence, on the magnitude of the quadratic phase distortion. Thus, satisfies the conditions
To simplify the subsequent analysis, we apply Eqs. (14) to a modified expression for . In accordance with Eqs. (11), the coordinates of the stationary point areFrom Eqs. (11), it follows thatChanging variables in Eqs. (16) by using Eqs. (15), we can obtainSolving the system of equations (17) for the second-order derivatives , , in terms of , , , and , and substituting the results into Eq. (13), finally simplifies the Hessian toWe take into account the fact that the MTF behavior, as seen from Eq. (10), is determined by , which, in turn, can be evaluated from Eq. (18):From Eq. (18), it is clear that Eqs. (14) are fulfilled if and are linear combinations of and ; in other words,where is the constant coefficient . Nonlinear solutions and of systems (11, 14) are beyond the scope of the present study. Here and in what follows, we limit ourselves to the analysis of phase masks described by polynomial expansions in x, y (see Appendix A for details).3. CUBIC PHASE MASK
Consider now the phase mask function chosen in the form of a third-order polynomial
where α is the magnitude of the phase deviation, a large real number, ; , , and , are the constants, and . In general, is an asymmetric function with respect to x and y. This, for example, is beneficial for imaging systems subject to asymmetric second-order aberrations. Substituting Eq. (21) into Eq. (7) results inFrom Eqs. (11, 22) it is clear that the function given by Eq. (21) ensures the required behavior defined by Eq. (20) for the stationary point coordinates and .For lower-order polynomials, for example , from Eqs. (11) it follows that there are no stationary points for the function . For higher-order polynomials, the stationary points are described by a system of two nonlinear equations (11), and the dependences of and on the parameters , 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 of the quadratic phase distortion only when the phase mask 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
and the system solution takes the form Employing Eqs. (24, 25), Eq. (19) simplifies toAs is seen from Eqs. (24, 25), the phase function has only one stationary point. Moreover, the denominators in Eqs. (24, 25, 26) are nonzero for all spatial frequencies only ifCondition (27) ensures that the quadratic form (in variables , ) represented by the denominator in Eqs. (24, 25, 26) is a sign definite quadratic form; i.e., it is always positive or negative (Sylvester criterion) [20].In the case of , from Eqs. (24, 25) it follows that the stationary point tends to the origin (, ), and, thus, the stationary point lies within the pupil area Ω. From Eq. (22), the Hessian matrix M, defined by Eq. (12), at becomes
The matrix signature is defined by the signs of the eigenvalues of the matrix M, which, in turn, are derived from the quadratic equationSo, the signature δ of M becomes a function of the spatial frequencies , and of the parameters , , , and of the phase mask.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 ) takes the form
As is seen from Eq. (30), the evaluated OTF is a well-defined and smooth function having no singularities at spatial frequencies . If , recalling Eq. (5), we obtain . Thus, the OTF of the imaging system with the phase mask given by Eq. (21), i.e., a generalized cubic mask, can be expressed by the following approximate equation (valid at all frequencies , ): Note that for practical purposes, the modulus of should be limited by unity at , .The precision of the OTF asymptotic representation, Eq. (30), at large α is determined by the residual term or . However, a higher precision can be achieved by evaluating additional terms of the order , where 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, , , and (ii) a cubic phase mask with no mixed terms, , .
In the first case, (i), Eq. (21) simplifies to , 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 , and at . One may conclude from Eq. (29) that the eigenvalues , of the Hessian matrix M have different signs at all , and, thus, . Analogously to Eq. (31), the OTF becomes
In the second case, (ii), the phase function is , which corresponds to the cubic phase mask originally suggested in [1]. From Eq. (29) , , and the signature of M evaluates to . Thus, Eq. (30) at and , gives
which for a square pupil (the pupil area in canonical pupil coordinates) coincides with the approximation obtained in [1], representing the defocused OTF of the imaging system with the cubic phase mask. At , the OTF approximation results inand, similarly, at the approximation of the OTF takes the formThe OTF given by Eq. (33) can be normalized to unity at , ; i.e., .4. DISCUSSION AND COMPARISON WITH DIRECT NUMERICAL SIMULATIONS
Asymptotic equation (10) presents the two-dimensional OTF of the incoherent fixed-focus imaging system subject to the quadratic phase distortion , given by Eq. (1), and having a phase-coding mask in its pupil plane. With the phase of the mask expressed by a third-order polynomial, Eq. (21), the OTF asymptotics at result in Eq. (30), which is invariant with respect to the magnitude of the phase distortion . 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 , , , and of the mask.
The derivation of Eqs. (10, 30) implies that the phase function 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 , and with the parameters of the mask satisfying Eq. (27), the OTF approximation is a smooth function with no singularities at spatial frequencies . In contrast to Eq. (30), Eq. (33) together with Eqs. (34, 35) exhibit singularities at and . 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., , evaluated with Eq. (31) and purely numerical calculations based on Fourier transformation of the pupil function with the cubic phase mask , defined by Eqs. (3, 21), respectively. The numerical calculations have been performed by using the OTF expressed via , for example [19],
where and denote the Fourier transformation pair in a symmetric form [18]. Equations (31, 36) have been simulated numerically in a C++ program. The program evaluates on a grid and utilizes the FFTW-3.2 library [22] to calculate the MTF in accordance with Eq. (36). For simplicity, we assume that the quadratic phase distortion is characterized only by the focusing error, i.e., .Figure 1 shows the MTF cross sections at calculated according to Eqs. (31, 36) for phase mask parameters and varying degree of defocus . Simulations have been carried out for three phase masks: plots 1–3, ; plots 4–6, ; plots 7–9, . In conformity with the stationary phase method, which requires , plots at show better agreement with the asymptotic curves than plots corresponding to . It can be also seen that at frequencies the phase functions with mixed terms and demonstrate better concordance with the approximations. One may conclude that a cubic phase filter with properly chosen and terms may eventually require a lower signal-to-noise premium in comparison with the phase mask function . 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 .
Often the phase mask function 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 of the quadratic phase distortion, defined by Eq. (1), only if the phase mask 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 specified by Eq. (1) the asymptotic representation of the OTF according to Eq. (10) becomes invariant with respect to the amplitude of only if the phase mask function is a third-order polynomial in x and y.
Invariance of on requires Eqs. (14) to be fulfilled. By recalling the Hessian definition, Eq. (13), Eqs. (14) can be rewritten in terms of , where and obey Eqs. (11),
or, after differentiation,Equations (A2) are fulfilled if and , 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 and, thus, second-order derivatives of F cannot be zeros, we obtainwhere are the coefficients. On the other hand, by substituting the Taylor expansions of and at into Eq. (7), it can be found thatBy equating Eqs. (A3, A4) and performing integration over x and y, we obtain the phase mask function in the form of a third-order polynomial (in x and y): being constants. The coefficients in Eq. (A3) can be expressed in terms of ,Finally, the linear expressions and corresponding to Eq. (20) can be obtained by substituting Eq. (A3) into Eq. (11).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.