An imaging design approach which is free of third-order astigmatism for one freeform optical surface and the image is presented in this paper. A set of differential equations is derived from generalized ray tracing. The solution of the above derived equations provides the anastigmatic freeform optical surface, the image surface and the object to image mapping. The obtained design can be used as a good starting point for optimization. As an example, a reflective freeform surface is designed for a single reflective Head Mounted Display (HMD). This example has a 3 mm pupil, 15mm eye clearance, 24-degree diagonal full field of view, and the final design yields an average MTF of 62.6% across 17 field points.
© 2014 Optical Society of America
Aspherical and freeform optical surfaces can provide more compact and lightweight imaging designs in most cases. Design methods based on multiparameter optimization depend on a good starting point to achieve a good local minimum of the merit function, and direct methods can provide a good starting point.
For rotational symmetrical aspherical designs with direct methods, Schmidt corrector plate is an early example of a single aspherical design to correct spherical aberration in a telescope , and Schwarzschild pioneered the design of two aspherical surfaces for full aplanatic system, i.e. they were free from spherical aberration and circular coma of all orders . Later, Wassermann and Wolf  generalized Schwarzschild approach presenting a method for the design of two aspherical surfaces as a solution of two simultaneous first-order ordinary differential equations. Standard integrating methods have been applied to obtain the numerical solution to get the aspherical profiles, although Willstrop and Lynden-Bell gave analytical solutions of the Schwarzschild design [4, 5]. The SMS 2D method is another method that has been used to directly design two rotational aspherical surfaces for imaging , which can be seen as a system of functional differential equations . It has been used, for instance, to design a telephoto for the SWIR band , and the SMS aspherical surfaces have been used as the starting point for finding an optimized solution [9, 10].
However, in the case of non-rotational symmetrical (freeform) surfaces for imaging, fewer direct methods have been proposed. A geometric algorithm has been proposed to directly design a single freeform surface for distortion correction  and a generalization of the Wassermann and Wolf differential equations has been suggested for two freeform surfaces . The SMS 3D method, which can be seen as a system of functional partial differential equations, has also been used for two freeform sharp-imaging (i.e., with no ray aberrations) of up to three object points, and also as a starting point for a posterior optimization .
In this paper we describe a novel direct method for designing a single aspherical or freeform surface which is free of astigmatism of third-order. In this method, only the object (or image) surface can be prescribed, and the image (object) surface is calculated, and the mapping between those surfaces is also obtained as a result of the design. The method is based on finding the differential equation on the optical surface deduced by equating the principal curvatures of the output wavefronts, which is done in Section 3 for both an aspherical rotational surface and a freeform surface. A brief description of getting the numerical solution is presented in Section 4. In Section 5, specific designs are shown, including the use of a freeform anastigmatic surface as a starting point for a subsequent optimization.
2. Statement of the problem
The method of single surface design presented here uses a biparametric set of rays prescribed at the object space defined with a wavefront as the input data, which will be refracted or reflected on the unknown refractive or reflective optical surface, and those output rays in the object space will be perpendicular to the output wavefront, which is unknown. We will refer to these two wavefronts as input and output “supporting” wavefronts, because for each ray of those wavefronts, we will consider another wavefront tangent to the supporting wavefront at that ray, which we will call “supported” wavefront, and whose principal curvatures and orientation of curvature lines do not coincide with that of the supporting wavefront. Therefore, we have a biparametric set of supported wavefronts tangent to one supporting wavefront. In our approach, we will consider a second order approximation of the supported wavefronts along the tangent point to the supported wavefront, and therefore, each supported wavefront will be fully characterized with the value of their principal curvatures and the orientation of the curvature lines relative to those of the supporting wavefront. This description is adequate, for instance, to model an optical system with a small aperture stop. A schematic drawing of the system is shown in Fig. 1 to clarify these definitions.
Considering that the input supporting wavefront and the second order approximations of the input supported wavefronts are given, the design problem is stated as follows: find the optical surface (refractive or reflective) such that the second order approximations of all the output supported wavefronts are spherical, so the design will be free from third-order astigmatism. The image surface will not be constrained, but calculated afterwards as the surface which the centers of curvature of the spherical output supported wavefronts lay on. These design conditions lead to a second-order partial differential equation in the general case, and to second-order ordinary differential equation in the rotational symmetrical case, as shown in the next section.
3. Differential equations for the optical surface profile
In this section, the analysis of both the freeform and the rotational symmetrical cases is for refractive surface. But it’s also valid for a reflecting optical surface.
3.1 The freeform case
The generalized ray tracing equations link the curvatures of the optical surfaces, the ones of the input and output wavefronts and the twist of the principal lines of curvature caused by refraction. Considering the plane of incidence and the plane perpendicular to the incident plane which contains the normal to the surface as well, the Coddington equations are the ones giving the propagation of the parameters of the supported wavefronts as [13, 14]:15]
We can describe the optical surface in spherical coordinates (r(θ, φ), θ, φ), where r is the radial distance, θ is the complement of the polar angle and measured from the plane to which the zenith direction is perpendicular, φ is the azimuthal angle. The first fundamental quantities of the optical surface are functions of r, θ, φ, ∂r/∂θ = rθ and ∂r/∂φ = rφ, while the second fundamental quantities are also functions of ∂2r/∂θ2 = rθθ, ∂2r/∂φ2 = rφφ and ∂r2/∂θ∂φ = rθφ, as follows:
We will restrict our analysis here to the particular case in which the input supporting wavefront is spherical centered at the origin, and the second order approximation of the input supported wavefronts are also spherical, so ρpi = ρqi = R(θ, φ) and 1/σi = 0.
Since the design condition is that the second order approximation of the output supported wavefronts must be spherical too, that is, ρpo = ρqo and 1/σo = 0, from Eq. (2) we deduce that 1/σs = 0, i.e., p and q are the principal directions for the optical surface. Thus the fundamental quantities F and M are 0 .
Extracting ρqo, ρpo from Eq. (1), the equation ρpo = ρqo gives as:
From Eq. (4) and Eq. (5) we can get the expressions of ρqs = E/L and ρps = G/N as functions of (r, θ, φ, rθ, rφ, rθθ, rφφ, rθφ), and therefore, Eq. (6) is a single equation linking (r, θ, φ, rθ, rφ, rθθ, rφφ, rθφ) with R, αi, αo, ni, and no. Since αi can be written as function of r, θ, φ, rθ, rφ and αo can be obtained with Snell’s law from αi also as function of r, θ, φ, rθ and rφ (ni, no are prescribed), Eq. (6) is a second-order Partial Differential Equation (PDE) with unknown function r(θ, φ).
A particularly simple expression is obtained if we restrict our analysis to second order approximation of the input supported wavefronts which are flat, i.e., R(θ, φ) = ∞. In this case, Eq. (6) is just
Direct computations of the ratio ρps/ρqs from the expressions obtained when solving Eq. (4) lead to ρps/ρqs = GL/EN, which can be expressed as a function of r, θ, φ, rθ, rφ, rθθ, rθφ rφφ using Eq. (5). On the other hand, αo can also be written as a function of r, θ, φ, rθ, rφ as shown in the Appendix A. Therefore, when we substitute E, G, L, N and αo(r, θ, φ, rθ, rφ), Eq. (7) is a second order PDE of the form:Eq. (1) which together with αo(θ, φ) will give the information to calculate the point on the image surface corresponding to the parameters (θ, φ).
3.2 The rotational symmetric case
In the particular case in which the function R(θ, φ) = R(θ), the input supporting wavefront and the input supported wavefronts are symmetric with respect to the axis θ = 0. If the boundary conditions of the PDE are also consistent with that symmetry, the solution will be found to be rotational symmetrical. In fact, such a solution can be obtained by a simpler second-order ordinary differential equation in r(θ), using that in this case the tangential and sagittal directions are the principal directions and therefore the expressions of ρps and ρqs are given by [16, 17]
Once r(θ) has been obtained, ρpo or ρqo can be calculated with Eq. (1) and with it can be calculated the image surface, which will be rotational symmetrical but not flat in general.
4. Numerical solution of PDE 
Following , the numerical solution to a hyperbolic PDE like Eq. (8) requires a boundary condition that can be set discretizing a curve with n points (, , …, ). Each one of them contains the values of r, θ, φ, rθ, rφ, rθθ, rφφ and rθφ that fulfill the PDE (Fig. 2 (a)). Through every point there will pass two distinct characteristic curves whose slopes are given as
The f-type characteristic of slope through the point will meet the g-type characteristic of slope through the point at the point with .straight line approximation. This gives the first approximation of and . by solving simultaneously
The first approximation of and can be obtained by solving simultaneously
The first approximation of is then given by
Now we solve the Eq. (12) for improved and with
Then the improved value of is found from Eq. (14). An iterative procedure is applied until all of r, θ, φ, rθ and rφ at the point have converged to (Fig. 2 (b)). Repeating the above procedure, we can obtain the numerical solution as a series of points (, , …, , …, , , …, , …, ). (Fig. 2 (c)).
5. Optical design example
Design examples of the rotational solutions have been presented previously , so only a freeform solution will be exemplified here.
The design is symmetrical about the plane φ = 0. Thus the intersection curves of the optical surface and the image surface with this plane is chosen as initial curves. The value of r, θ, rθ and rθθ can be obtained from the design presented previously . From the symmetry, we can deduce that φ = 0, rφ = 0 and rθφ = 0. rφφ can be obtained by substituting the above values into the PDE. Thus they fulfill the condition that all the output supported wavefronts converge to the point on the initial curve of the image surface (Fig. 3). From the initial curve of the optical surface, other points of the optical surface have been obtained by solving the PDE applying the numerical method described in Section 4. After that the optical surface is represented using a Forbes polynomial fitted to these points.
The example of a freeform anastigmatic reflector design is shown in Fig. 4. The whole system has 15.5 mm eye clearance, 24 degree diagonal full field of view (9.56° x semi-field and 7.2° y semi-field), and a 14° mirror tilt angle. The distance from the pupil to the vertex of the mirror is 16.0mm. Distance from the vertex of the mirror to the image plane is kept around 13mm, The image plane has a rectangular aperture with a size of 4.8mm by 3.6mm in the x and y dimensions, respectively. The effective focal length is 13.2mm and the f-number is 4.4.
After obtaining the numerical solution with the prescribed boundary condition, we have calculated the two surfaces formed by the centers of principal curvatures of the output supported wavefronts separately. Both surfaces should be coincident, so the differences between them is a measure of the calculation errors. The max deviation is 0.01% of the average radius, i.e., it is essentially anastigmatic.
This design method has not enough degrees of freedom to prescribe a flat image surface, which, of course, is much more interesting than a curved one for practical applications. Nevertheless, the boundary conditions can be chosen to get an “almost” flat image surface. For this purpose we have chosen a straight line as the initial curve of the image surface. A subsequent optimization of the reflector surface has been done in CodeV evaluating the image quality at the plane which best approximates the curved image surface. The result of this optimized design is shown in Fig. 5. The Modulation Transfer Function has been evaluated for 17 field points and the result is shown in Fig. 5(left). The distortion grid is plotted in Fig. 5 (right).
Cakmakci et al  have proposed a combination of Gaussians as a sum of basis functions proposed for the representation and optimization of the optical surface profile in a similar configuration with 60.5% of average MTF value at the frequency of 23 cycles/mm and 3.6% of max distortion. The design example presented in this paper has achieved a better average MTF value of 62.6% at the frequency of 23 cycles/mm. The max distortion is 9.65%.
6. Discussion and conclusion
Since the perfect imaging requires that the tangential foci and sagittal foci are the same, we have applied this differential equation approach in imaging designs by considering the tangential and sagittal ray propagation simultaneously. This means that the freedom is removed from controlling the mapping from the object to image, resulting that the mapping from the object to image as well as the shape of the image surface is obtained after the design of the optical surface. But instead, it is used to control the sagittal ray propagation, providing the significant gain in image quality evaluations such as MTF value.
Based on this approach, a single surface anastigmatic design example has been developed. As a good initial design for optimization, a single mirror HMD design with an average MTF value of 62.6% at the frequency of 23 cycs/mm has been achieved.
Authors thank the European Commission (SMETHODS: FP7-ICT-2009-7 Grant Agreement No. 288526, NGCPV: FP7-ENERGY.2011.1.1 Grant Agreement No. 283798), the Spanish Ministries (ENGINEERING METAMATERIALS: CSD2008-00066, SEM: TSI-020302-2010-65 SUPERRESOLUCION: TEC2011-24019, SIGMAMODULOS: IPT-2011-1441-920000, PMEL: IPT-2011-1212-920000), UPM (Q090935C59) and the academic license for CodeV from Synopsys for the support given to the research activity of the UPM-Optics Engineering Group, making the present work possible.
References and links
1. B. Schmidt, “Ein lichtstarkes komafreies Spiegelsystem,” Central-Zeitung für Optik und Mechanik 52(2), 25–26 (1931).
2. K. Schwarzschild, “Astronomische Mitteilungen der Königlichen Sternwarte zu Göttingen 10, 3 (1905),” Reprinted: Selected Papers on Astronomical Optics. SPIE Milestone Ser. 73(3), (1993).
3. G. D. Wasserman and E. Wolf, “On the theory of aplanatic aspheric systems,” Proc. Phys. Soc. B 62(1), 2–8 (1949). [CrossRef]
4. D. Lynden-Bell, “Exact optics: a unification of optical telescope design,” Mon. Not. R. Astron. Soc. 334(4), 787–796 (2002). [CrossRef]
5. R. V. Willstrop and D. Lynden-Bell, “Exact optics – II. Exploration of designs on- and off-axis,” Mon. Not. R. Astron. Soc. 342(1), 33–49 (2003). [CrossRef]
6. R. Winston, J. C. Miñano, and P. Benítez, Nonimaging Optics (Academic, 2005), Chap. 9.
7. F. Duerr, P. Benítez, J. C. Miñano, Y. Meuret, and H. Thienpont, “Analytic free-form lens design in 3D: coupling three ray sets using two lens surfaces,” Opt. Express 20(10), 10839–10846 (2012). [CrossRef] [PubMed]
8. J. M. Infante, “Optical Systems Design using SMS method and optimizations,” Ph.D. Thesis, Universidad Politécnica de Madrid (2013).
11. J. Hou, H. Li, Z. Zheng, and X. Liu, “Distortion correction for imaging on non-planar surface using freeform lens,” Opt. Commun. 285(6), 986–991 (2012). [CrossRef]
12. D. Cheng, Y. Wang, and H. Hua, “Free form optical system design with differential equations,” Proc. SPIE 7849, 78490Q (2010). [CrossRef]
13. O. N. Stavroudis, The Mathematics of Geometrical and Physical Optics: The k-function and its Ramifications (Wiley-VCH, Berlin, 2006).
14. C. Wang and D. L. Shealy, “Multimirror anastigmat design,” Proc. SPIE 2000, 24–33 (1993). [CrossRef]
15. D. J. Struik, Lectures on Classical Differential Geometry: Second Edition (Dover Books on Mathematics, 1988).
16. A. Estrada-Molinam and R. Díaz-Uribe, “Tangential and sagittal curvature from the normals computed by the null screen method in corneal topography,” Proc. SPIE 8011, 80119J (2011). [CrossRef]
17. E. Weisstein, “Radius of Curvature.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/RadiusofCurvature.html
18. E. H. Twizell, “The numerical solution of second order hyperbolic partial differential equations with unequally spaced initial conditions,” Comput. J. 18(3), 252–257 (1975). [CrossRef]
19. J. Liu, J. C. Miñano, P. Benítez, and L. Wang, “Single optical surface imaging designs with unconstrained object to image mapping,” Proc. SPIE 8550, 855011 (2012). [CrossRef]
20. O. Cakmakci, B. Moore, H. Foroosh, and J. P. Rolland, “Optimal local shape description for rotationally non-symmetric optical surface design and analysis,” Opt. Express 16(3), 1583–1589 (2008). [CrossRef] [PubMed]