Orthogonal basis with a conicoid first mode for shape specification of optical surfaces

A rigorous and powerful theoretical framework is proposed to obtain systems of orthogonal functions (or shape modes) to represent optical surfaces. The method is general so it can be applied to different initial shapes and different polynomials. Here we present results for surfaces with circular apertures when the first basis function (mode) is a conicoid. The system for aspheres with rotational symmetry is obtained applying an appropriate change of variables to Legendre polynomials, whereas the system for general freeform case is obtained applying a similar procedure to spherical harmonics. Numerical comparisons with standard systems, such as Forbes and Zernike polynomials, are performed and discussed. © 2016 Optical Society of America OCIS codes:(220.1250) Aspherics; (220.4830) Optical systems design; (220.4840) Optical testing; (220.4610) Optical fabrication; (000.4430) Numerical approximation and analysis. References and links 1. R. N. Wilson,Reflecting Telescope Optics I , (Springer, 2004). 2. Y. Ding, X. Liu, Z.-R. Zheng, and P.-F. Gu, “Freeform LED lens for uniform illumination,” Opt. Express 16, 12958–12966 (2008). 3. M. Victoria, C. Domı́nguez, I. Antón, and G. Sala, “Comparative analysis of different secondary optical elements for aspheric primary lenses,” Opt. Express 17, 6487–6492 (2009). 4. Y. Wu, J. Xi, M. J. Cobb, and X. Li, “Scanning fiber-optic nonlinear endomicroscopy with miniature aspherical compound lens and multimode fiber collector,” Opt. Lett. 34, 953–955 (2009). 5. D. J. Meister and S. W. Fisher, “Progress in the spectacle correction of presbyopia. Part 2: Modern progressive lens technologies,” Clin Exp Optom. 91, 251–264 (2008). 6. R. Navarro, “Adaptive model of the aging emmetropic eye and its changes with accommodation,” J. Vis. 14(21), 1–17 (2014). 7. F. Chen, G. M. Brown, M. Song, “Overview of three-dimensional shape measurement using optical methods,” Opt. Engineering39, 10–22 (2000). 8. D. Michaelis, P. Schreiber, and A. Bräuer, “Cartesian oval representation of freeform optics in illumination systems,” Opt. Lett.36, 918–920 (2011). 9. H. Ries and J. Muschaweck, “Tailored freeform optical surfaces,” J. Opt. Soc. Am. A 19, 590–595 (2002). 10. G. W. Forbes, “Shape specification for axially symmetric optical surfaces,” Opt. Express 15(8), 5218–5226 (2007). 11. A. B. Bhatia, E. Wolf, and M. Born, “On the circle polynomials of Zernike and related orthogonal sets,” Proc. Camb. Philos. Soc. 50(1), 40–48 (1954). 12. G. W. Forbes, “Robust, efficient computational methods for axially symmetric optical aspheres,” Opt. Express 18(19), 19700–19712 (2010). #249995 Received 16 Sep 2015; revised 18 Nov 2015; accepted 1 Dec 2015; published 3 Mar 2016 (C) 2016 OSA 7 Mar 2016 | Vol. 24, No. 5 | DOI:10.1364/OE.24.005448 | OPTICS EXPRESS 5448 13. G. W. Forbes, “Manufacturability estimates for optical aspheres,” Opt. Express 19(10), 9923–9941 (2011). 14. G. W. Forbes, “Characterizing the shape of freeform optics,” Opt. Express 20(3), 2483–2499 (2012). 15. G. W. Forbes, “Fitting freeform shapes with orthogonal bases,” Opt. Express 21(16), 19061–19081 (2013). 16. R. Navarro, J. L. López, J. A. Dı́az, and E. Pérez Sinusı́a, “Generalization of Zernike polynomials for regular portions of circles and ellipses,” Opt. Express 22, 21263–21279 (2014). 17. C. Ferreira, J. L. López, R. Navarro , and E. Pérez Sinusı́a, “Zernike-like systems in polygons and polygonal facets,” Appl. Opt.54, 6575–6583 (2015). 18. E. W. Weisstein, “Legendre Polynomial”, From MathWorld. A Wolfram Web Resource. http://mathworld.wolfram.com/LegendrePolynomial.html . 19. E. W. Weisstein, “Spherical Harmonic”, From MathWorld. A Wolfram Web Resource. http://mathworld.wolfram.com/SphericalHarmonic.html . 20. E. W. Weisstein, “Associated Legendre Polynomial”, From MathWorld. A Wolfram Web Resource. http://mathworld.wolfram.com/AssociatedLegendrePolyno mial.html.


Introduction
The number and relevance of applications of aspheric and freeform optics is continuously increasing, ranging from astronomy [1], industry [2], solar energy [3], biomedical optics [4], or physiological optics [5], among others. The high complexity of optical surfaces found in biological systems such as the human eye [6], or the new advances in fabrication and testing of freeform surfaces [7], are demanding precise, robust and efficient methods of specifying these surfaces. Ideally, the shape specification should be physically meaningful and invariant for the different stages of design, fabrication, testing or application. In optics, it is common to represent the surface sag z as a function of the coordinates z = f (x, y). Often the configuration is nearly rotationally symmetric and hence it is better to work in cylindrical coordinates z = f (r, θ ). The variety of systems of representations range from the pure sampling grid of points, localized splines, or global or modal representations given by combinations of functions such as spheres, conicoids, monomials, polynomials, etc. Specific methods of representation tailored for specific applications were also proposed, such as generalized Cartesian ovals [8] or solutions of specific differential equations [9], among others. Due to the high relevance of spheres, the most widely used characterizations of optical surfaces is the sum of a sphere (or conicoid C) plus an aspherical part A usually given as a linear combination of terms z = C + A. In what follows we will talk of conicoids C which include the sphere as a particular case.
Typically the terms specifying A are either monomials or polynomials. Historically, monomials were used first due to their apparent simplicity, but as Forbes pointed out referring to monomials [10], "the most widely used characterization of surface shape is numerically deficient", mainly due to their lack of orthogonality. Orthogonal systems of polynomials, such as Zernike polynomials [11], Forbes polynomials [10,[12][13][14][15], etc., permit to overcome a series of key issues ranging from numerical instabilities to effective tolerance specifications. Furthermore, basic linear algebra tells us that the two crucial properties of a good system of representation (sets of basis functions) are orthogonality and completeness. An additional, but less crucial property is normality (norm unity for all basis functions). Orthogonality, that implies the linear independence between the basis functions, implies also good numerical behavior, avoids redundancy and ensures uniqueness of the representation, among other highly important properties. In this context completeness is even more crucial as it means that the system can represent all possible surface shapes, that is to insure that we have a real freeform system.
In this sense the most widely used methods of optical surface representation in the form of z = C + A are essentially non orthogonal, even when they use Zernike polynomials, Forbes polynomials or orthogonal systems to represent the departure from the sphere (or conicoid), simply because the conicoid itself is not orthogonal to A. Here our goal was to solve this problem to obtain a system of representation in which C is orthogonal to A (of course the basis functions of A are orthogonal as well). This implies that C is one of the basis functions of the system.
To this end, we divided the main goal into the following specific objectives: (1) to develop a general theoretical framework to obtain this type of systems, initially restricted to rotationally symmetric surfaces z = f (r); (2) to obtain and implement a complete orthogonal system, in which the first basis function is a conicoid; (3) to generalize the above results to non-rotationally symmetric surfaces which is essential to obtain true freeform systems; and (4) to test numerically our new systems and perform direct comparisons with previous methods. The next Sections are organized accordingly.

Basis for rotationally symmetric surfaces
In this section we introduce the general framework of our theory, restricted to rotationally symmetric surfaces defined by an equation of the form z = f (r), r ∈ [0, 1], where r and z are cylindrical coordinates. We design an orthogonal system for L 2 ν [0, 1] with measure dν = rdr, in which the first element of the system is a specified conicoid C that, in the following, we denote by q 0 (r) for convenience. The remaining elements of the system, the elements of the set A, are denoted by q n (r), n = 1, 2, 3, . . . The functions q n (r), n = 1, 2, 3, . . ., are constructed using three essential ingredients: (i) an arbitrary orthonormal system {p n (x)} n=0,1,2,... , with respect to a certain measure dµ = ρ(x)dx in an interval [c, d], (ii) the selected conicoid q 0 (r) and (iii) a convenient change of variable x = ϕ(r), ϕ : A similar method based on a change of variables was successfully applied before to obtain orthogonal Zernike-like sytems on noncircular apertures [16], polygons and polygonal facets [17]. The resulting orthogonal system consists of functions {q 0 (r), q 1 (r), q 2 (r), . . .} defined in the interval [0, 1] that are orthogonal with respect to the measure dν = rdr. Moreover, the functions q 1 (r), q 2 (r), . . . have also norm unity. Rougly speaking, the idea is the following: we use ϕ(r) to replace the first element p 0 of the system {p n (x)} n=0,1,2,... by q 0 (r). To preserve the orthogonality of the new system, we must choose ϕ(r) appropriately. In the remaining of this section we develop this idea and give an important example.
We observe that, when we define q n (r) := C n p n (ϕ(r))w(r), w(r) := ρ(ϕ(r))ϕ ′ (r) with C n = 1 for n = 1, 2, 3, . . ., and C 0 arbitrary at this moment; we find that the functions q n (r), n = 1, 2, 3, . . . are orthonormal in [0, 1] with respect to the measure dν = rdr. It is worth to note that the variable r in the denominator of w(r) is not dangerous as the numerator behaves as r when r → 0, as we see in the first line of (4). To complete our basis, we need to introduce a first element q 0 (r), the arbitrary conicoid, assuring that it is orhogonal to all the elements q n (r), n = 1, 2, 3, . . . This is achieved by choosing the function ϕ(r) as the unique solution of the boundary value problem: Since the first-order differential equation in Eq. (4) can be directly solved by integration, the unique solution x = ϕ(r) of Eq. (4) is implicitly defined by the equation when we take It is obvious that the left hand side of Eq. (5) is an increasing function of x and the right hand side of Eq. (5) is an increasing function of r. Then, ϕ(r) is a monotonic function with ϕ ′ (r) > 0. Thus, we have that the set {q n (r)} n=0,1,2,... is a quasi-orthonormal basis of L 2 ν [0, 1] (orthonormal except for the fact that ||q 0 || 2 r = C 2 0 ) with q n (r) = q 0 (r) C 0 p 0 C n p n (ϕ(r)), n = 0, 1, 2, . . . Moreover, that belongs to L 2 µ [c, d]: with

Solution for normalized Legendre polynomials
The normalized Legendre polynomials [18] p n (x) = 2n + 1 2 are an orthonormal basis of L 2 [−1, 1] with respect to the weight function ρ(x) = 1. We take the following conicoid: with s = ±1. If s = 1 and b > 0, q 0 is the portion of an ellipse with semiaxes a and b in the upper half-plane between the angles arctan(b √ a 2 − L 2 /a) and π/2, and with 0 < L < a. If s = 1 and b < 0, q 0 is the portion of an ellipse with semiaxes a and −b in the lower halfplane between the angles −π/2 and − arctan(b √ a 2 − L 2 /a), and with 0 < L < a. If s = −1 and b > 0, q 0 is the portion of a hyperbola with semiaxes a and b in the upper half-plane between the angles arctan(b √ a 2 + L 2 /a) and π/2, and if s = −1 and b < 0, q 0 is the portion of a hyperbola with semiaxes a and −b in the lower half-plane between the angles −π/2 and − arctan(b √ a 2 + L 2 /a). In any case, from Eq. (5), We obtain the value of C 0 from Eq. (6), Thus, Therefore, we have that the set is a quasi-orthonormal basis of L 2 ν [0, 1] with dν = rdr and any function F(r) ∈ L 2 ν [0, 1] can be written as with c n , C 0 and q n (r) given in Eqs. (11), (16) and (18) respectively. In Table 1 we can find the first five functions q 0 (r), q 1 (r), . . . , q 4 (r). The graphs of these functions (up to n = 4) are illustrated in Fig. 1 for a elliptical cap, r ∈ [0, 1], s = 1, b = 1 and a = L = 3/4.

Observation 1.
In optics, it is standard to express rotationally symmetric surfaces in terms of their deviation from the sagittal representation. If we consider as first term (18)) obtained from the normalized Legendre polynomials and q 0 ( with 0 < L 2 ≤ 1/(εc 2 ) if ε > 0, and as initial orthonormal basis the normalized Legendre polynomials given in Eq. (13), we obtain a new quasi-orthonormal basis with and that is, This is a significantly less compact expression as compared to Eq. (18).
The theory developed in this section only applies to rotationally symmetric surfaces specified by an equation of the form z = f (r), with r ∈ [0, 1]. In the following section we generalize the theory to arbitrary surfaces defined in cylindrical coordinates (r, θ , z), specified by an equation of the form z = f (r, θ ).

Basis for freeform surfaces
In this section we formulate a more general theory to approximate arbitrary surfaces z = f (r, θ ) defined over the unit disk (r cos θ , r sin θ ) ∈ D. We design an orthogonal system for L 2 ν (D) with measure dν = rdrdθ in which the first element of the system is a specified rotationally symmetric surface q 0 0 (r). The remaining elements of the system q m n (r, θ ), n = 1, 2, 3, . . ., are constructed using again the three essential ingredients used in the previous section: (i) an arbitrary orthogo- with C m n = 1 for (n, m) = (0, 0) and C 0 0 arbitrary at this moment, we find that the functions q m n (r, θ ), n, m(n) = 0, 1, 2, . . ., are orthonormal in D with respect to the measure dν = rdrdθ . To complete our basis, we need to introduce a first element q 0 0 (r) assuring that it is orthogonal to all the other elements q m n (r, θ ), n, m(n) = 0, 1, 2, . . ., (n, m) = (0, 0). This is achieved by choosing the function ϕ(r) the unique solution of the boundary value problem Eq. (4) given by Eq. (5).

Implementation and examples
For the numerical implementation and testing of the new basis proposed in Sections 2 and 3, we consider a Gaussian surface, as it is expected to require higher order expansions to obtain reasonable accuracies. In the first test, we use a rotationally symmetric Gaussian surface to compare with Forbes polynomials for the 1D case. The second test considers the same symmetric Gaussian surface to establish a comparison with our new basis in 2D, and the basis in 1D and Zernike polynomials. Finally, in the third test we use a non-symmetric elliptic Gaussian surface to compare our 2D system with Forbes and Zernike polynomials. In all the examples, we use the least square approximation for computing the coefficients.

Example 1.
In order to check the accuracy of the approximation supplied by the basis {q n (r)} n=0,1,2,... given in Eq. (18) and compare with Forbes' approximations, we consider the following Gaussian surface of revolution (see Fig. 3), For Forbes' approximation, we implement both, Q con m and Q bfs m polynomials [10,12]. On the one hand, Forbes' approximation using Q con m polynomials [10] reads where c and ε are, at this moment, free parameters, u is the normalized radial coordinate given by u = r/r max , r max = L is the aperture size, and Q con m (x) = P (0,4) m (2x − 1) with P (α,β ) m (x) the Jacobi polynomials of parameters (α, β ) = (0, 4). Following Forbes' algorithm, the values of c and ε in the first term on the right hand side of Eq. (36) can be approximated by c ≃ −0.398942 and ε ≃ −8.646945. Then, the first term in Eq. (36) is a hyperbola. On the other hand, we consider Q bfs m polynomials [10,12]. In this case where c bfs is the curvature of the best-fit sphere, r max = L and u = r/r max . Polynomials Q bfs m (x) can be generated using a non-standard recurrence relation [12,15] that involves a set of orthogonal polynomials. For this example, the best-fit curvature c bfs can be approximated by −0.306394.
Finally, we consider our expansion defined by Eq. (19), with {q n (r)} n=0,1,... the new quasiorthonormal basis given in Eq. (18), C 0 = 1 2 [2b 2 − s b 2 a 2 L 2 ] 1/2 and c n defined in Eq. (11). As the first order approximation is given by Eq. (14), in order to get a first close-fitting approximation, we must translate the original function, that is, f (r) + b, in such a way that both functions are coincident at r = 0. For this function, we consider s = 1, a = b = 1/c bfs in Eq. (14). The fit error obtained with the three approximations Eq. (19) (proposed here), Eq. (36) (Forbes Q con m polynomials) and Eq. (37) (Forbes Q bfs m polynomials) for the same number of terms (N = 7, M = 6 and M = 6 respectively, that is, 8 terms for the three cases) and for L = r max = 1, are compared in Fig. 4.
In Fig. 5 we compare the relative error provided by equations Eq. (19) (new basis proposed in Section 2) and Eqs. (36) and (37) (Forbes polynomials) in the L 2 −norm for the same number of terms. Our basis functions provide lower errors than Forbes Q con m polynomials, but not with respect to Forbes Q bfs m polynomials.

Example 2.
In order to check the accuracy of the approximation supplied by the basis {q m n (r, θ )} n,m(n)=0,1,2,... for the freeform case given in Eq. (34), we implement the same example as in the 1D case, but now we consider the approximation problem in 2D, We approximate this function by using the new quasi-orthonomal basis {q m n (r, θ )} n,m(n)=0,1,... given in Eq. (34) but choosing q 0 0 (r) in Eq. (14) the sphere with s = 1, a = b = 1/c bfs and L = 1, q 0 0 (r) = (c bfs ) −2 − r 2 : with C 0 = π 2 2c −2 bfs + 1 1/2 and c m n defined in Eq. (31). Here we compare with the approximation obtained with Zernike polynomials (ZPs) [11] f (r, θ ) ≃ and with the approximation obtained with Forbes polynomials (37). The results are similar to the ones obtained in Fig. 5 for the 1D case in comparison with Forbes polynomials. Our method compares favorably with ZPs for 2D.
We consider now the non-symmetric two-dimensional Gaussian type surface (see Fig. 6), In Fig. 7   a ≃ 0.653754 and b ≃ 1.02225. Moreover, using the idea given by Forbes, we also approximate the function using the new quasi-orthonomal basis in Eq. (34) where now q 0 0 (r) is the best-fit sphere, and then s = 1, and a = b = (0.261082) −1 . We also compare with ZPs Eq. (40). The comparison is given using the number of terms (calculated coefficients) where we take the best RMSE given by Forbes approximation (42) for different values of N and M = 0 up to 2N + M = 10. We also compare with the ZPs approximation. Our basis functions provide lower errors in comparison with Forbes approximation, and the difference decreases with the number of terms considered. In comparison with ZPs approximation, our basis funcions provide lower errors up to 10 terms, and ZPs provides lower errors from 15 terms, and also the difference decreases with the number of terms considered.

Discussion and conclusions
In summary, we developed a general, rigorous and powerful framework to obtain orthogonal systems to represent freeform optical surfaces. The method consists of first selecting the desired first basis function q 0 , a sphere or conicoid in our case, and then choosing an appropriate orthogonal system on the desired domain or support. For the case of rotational symmetry we considered the interval [0, 1] and the unit disk in general. For these 1D and 2D domains we selected Legendre polynomials and spherical harmonics respectively, as the initial systems. Then the method consists of finding a change of variables ϕ which transforms these systems {p n } n=0,1,2,... into another orthogonal system {q n } n=0,1,2,... in which the first function q 0 is the conicoid that we chose. Note that this theoretical framework is powerful enough to obtain orthogonal systems for different expressions of the conicoids: canonic expression of Eq. (14), or the sag of Eq. (20) commonly used in optical design and testing. We also included the semidiameter L of the surface explicitly, which permits to avoid the need of normalizing the radial coordinate r, that is necessary when using ZPs or similar systems. For the implementation and examples we used the canonical expression for the sake of simplicity, since it allows more compact expressions as compared to the sag equation. Compared to other systems used in optics, this work is more general in different aspects. Probably, ZPs form the most widely used orthogonal system in optics. It is the standard basis for representing wave aberration (or optical path difference), and it is also used for representing aspherical and freeform surfaces. This system includes the paraboloid of zero mean (defocus term Z 0 2 = 2r 2 − 1, r ∈ [0, 1]), but it requires to apply a higher order expansion to approximate spheres or other conicoids. In fact, it is common to use ZPs, Z m n to represent only the aspherical terms in the form z = C + A where C is the conicoid and A = ∑ n,m a m n Z m n ; a m n are the Zernike expansion coefficients. The approach developed by Forbes [10] (actually he proposed two different approaches for mild and strong aspheres respectively) is similar in the sense that it splits the surface sag into the same conic and aspherical parts z = C + A, but applying a smarter ad hoc expansion A = r 4 ∑ n f n Q n (r 2 ) where Q n are orthogonal. Among other advantages, extracting the common factor r 4 permits to reduce the order (and number) of Forbes polynomials Q n needed to approximate A. Nevertheless, the main difference of our approach is that the conicoid is the first element of the system and hence it is orthogonal with the rest. In this sense we propose a fully orthogonal system as opposed to these partially orthogonal systems applied only to the aspherical part A.
For the numerical implementation and testing we have chosen a Gaussian as it is expected to require higher order expansions (high values of n) to obtain reasonable accuracies. Our results with 1D (Legendre polynomials) and 2D (spherical harmonics) implementations provide results which compare favorably with Forbes polynomials and ZPs in the general freeform two-dimensional case, whereas Forbes Q bfs polynomials show a better performance in the one-dimensional rotationally-symmetric case. We want to remark that the cases and implementations presented here are particular examples of a much more general theoretical framework. In fact, it may be possible to consider a wide variety and types of initial surface q 0 as well as to use different types of polynomials, etc. The main restriction to obtain systems with analytical expressions is that the integral equations, Eqs. (4) and (5), that we have to solve to obtain the change of variables, must have an analytical solution. Nevertheless, we hope that the examples of conicoid-based orthogonal systems presented here are useful and general enough for designing, manufacturing and testing freeform optical surfaces.