Fitting discrete aspherical surface sag data using orthonormal polynomials

David Hilbig; Ufuk Ceyhan; Thomas Henning; Friedrich Fleischmann; Dietmar Knipp

doi:10.1364/OE.23.022404

1. Introduction

Nowadays, optical components with at least one aspherical surface are in widespread use. Most of these surfaces are rotationally symmetric with weak to strong deviation from a basis spherical shape. Modern manufacturing technologies are able to produce aspherical surfaces with sub-micrometer precision. Reliable testing of such components during fabrication is mandatory to assure the demanded quality. The outcome of a typical surface measurement is a series of discrete data points with inherent measurement errors. Analysis of the test result in order to obtain a meaningful characterization of the device-under-test usually involves finding the best-fit of a model description of the aspherical surface to the discrete data set. This is especially true for experimental ray tracing [1], where the surface shape can be reconstructed from ray slope data taken from a transmission measurement by means of a minimization process whose core is the used model for the surface description. Such a surface could be described using various common polynomials, as Zernike polynomials [2], Chebyshev polynomials [3] or splines [4]. At this point, the choice of the mathematical description can greatly improve this process. Forsythe [5] presented a way to greatly reduce the effort to solve the problem of the minimization process using the inner product on orthogonal polynomials which can be even further simplified by using orthonormalized polynomials. Section 2 describes the minimization problem, the simplified solution along with the common method using direct matrix inversion.

The rotational symmetry around a common axis in most aspherical surfaces lend itself to the application of cylindrical coordinates with an independence of the azimuthal angle φ. Describing the surface's shape with respect to its axial tangent plane z = f(ρ) leads to a sagittal representation which is commonly referred to as the surface sag. The most commonly used representation of this kind is the standard description for rotationally symmetric aspherical surfaces according to ISO 101110-Part12 [6], where the total asphere is decomposed into a base conic component and an additional power series for the higher order deviations from the conic.

z (ρ) = \frac{C ρ^{2}}{1 + \sqrt{1 - C^{2} ρ^{2} (1 + κ)}} + \sum_{m = 0}^{M} a_{m} ρ^{2 m + 4}

Here, C = 1/R is the paraxial surface curvature and κ is the conic constant defining the type of conic section. However, the power series used in this representation was reported by Forbes and Brophy [7] to be highly inefficient, numerical unstable and prone to round-off errors. One mayor reason for this is, as will be shown later, that this polynomial set is not orthogonal. Hence, the individual terms influence each other leading to a strong cancellation between the terms where the total number of polynomial orders will have a profound influence on each coefficient of the polynomial set in a fitting situation. An alternative representation that should overcome the drawbacks of this description was presented by Forbes along with a presentation of its enhanced performance in design and testing [8]. He introduced two new surface descriptions designated as Q-polynomials where the main focus was set to provide a representation containing an orthogonal set of polynomials. One of them was specially created with interferometrical testing in mind, where the description allows to constrain the maximal slope of the normal departure from a best-fit sphere over the coefficients of the polynomial expansion. The set was later expanded to yield a gradient-orthogonal basis for the description of freeform surfaces [9]. Its chosen orthogonality plays well to interferometric testing. In addition, slope-orthogonal polynomials were reported to excel in the minimization of typical optical performance parameters as the RMS spot size of a focusing system in a design situation [10]. However, this work will focus on the generally less regarded sag orthogonal basis that was referred to as the representation more suitable for general-purpose aspheres. Their close relation to the standard ISO-description and especially their orthogonality in the surface makes them better suited for the aforementioned application in aspherical surface reconstruction by experimental ray tracing. Another example are tactile systems which will not benefit from a description orthogonal in the surface slopes. Using finite differences to compute these from the discrete surface data will introduce round-off errors that overcome the advantages of an orthogonal set. A description of the Q_con polynomial set as well as an explanation of an efficient computation method to evaluate these to arbitrary high orders will be given in section 3. Similar to case of orthogonal Zernike polynomials as stated by Malacara and DeVore [11] as well as Mahajan [12], such a polynomial set, though orthogonal in the continuous sense, is not orthogonal anymore in case of discrete data sets. This is especially true for data sets with smaller numbers of data points which is a common situation found with surface testing based on Shack-Hartmann or the experimental ray tracing. Therefore, a Gram Schmidt based orthonormalization process [13, 14] is proposed in order to obtain a perfect orthogonal polynomial set for fitting to discrete surface data sets. As was shown lately by Ye et al., a similar approach could be successfully applied to the analysis of wavefront data over generally shaped apertures using Zernike polynomials [15].

2. Least-squares minimization with polynomials

Obtaining the best-fit of a model function S to a given set of measured data S_n with N data points usually involves a minimization process where the difference between both is minimized. Using polynomials the function S can be expanded by a polynomial P of degree L as

S = \sum_{l = 1}^{L} b_{l} P_{l} .

The best-fit to S_n is found by choosing the coefficients b_l to minimize

E^{2} (b_{1}, b_{2}, ..., b_{L}) = \sum_{n = 1}^{N} {[S_{n} - \sum_{l = 1}^{L} b_{l} P_{l} (x_{n}, y_{n})]}^{2} = minimum .

Here, E² is the error that needs to be minimized and represents a differentiable quadratic function of the coefficients b_l. Limiting values for b_i in an interval of real numbers where this function is continuous assures the existence of a minima according to the extreme value theorem. With E as a differentiable function, a minimum will be found if the partial derivatives of the function with respect to each coefficient are zero,

\frac{\partial E^{2}}{\partial b_{l}} = 0.

which means for Eq. (3)

\frac{\partial E^{2}}{\partial b_{l}} = 2 \sum_{n = 1}^{N} [S_{n} - \sum_{m = 1}^{M} b_{m} P_{m} (x_{n}, y_{n})] (- P_{l} (x_{n}, y_{n})) = 0

for all l = 1, …, L and M = L. Resolving the brackets and cancelling 2 leads to

\frac{\partial E^{2}}{\partial b_{l}} = \sum_{m = 1}^{M} b_{m} \sum_{n = 1}^{N} P_{m} (x_{n}, y_{n}) P_{l} (x_{n}, y_{n}) - \sum_{n = 1}^{N} S_{n} P_{l} (x_{n}, y_{n}) = 0,

so that at each minimum, one receives

\sum_{m = 1}^{M} b_{m} [\sum_{n = 1}^{N} P_{m} (x_{n}, y_{n}) P_{l} (x_{n}, y_{n})] = \sum_{n = 1}^{N} S_{n} P_{l} (x_{n}, y_{n}) .

With

c_{l} = \sum_{n = 1}^{N} S_{n} P_{l} (x_{n}, y_{n}) .

and the square coefficient matrix

G_{m, l} = \sum_{n = 1}^{N} P_{m} (x_{n}, y_{n}) P_{l} (x_{n}, y_{n}) .

Equation (7) changes to

\sum_{m = 1}^{M} b_{m} G_{m, l} = c_{l},

which represent a system of M equations with M unknowns that must be satisfied by the coefficients b_m at any minimum. These are the normal equations of the least-squares data fitting problem where a unique set of solutions exist if the determinant of G is unequal to 0. A classical numerical technique to solve this would be the Gaussian elimination method whose description can be found in standard textbooks about linear algebra in statistics [16]. However, a more elegant method can be found using matrix representations.

Equation (2) can as well be represented more conveniently in a matrix form as

S = b P .

Expanding this with the transpose of the polynomials, one gains the matrix equivalent of Eq. (7),

P^{T} P b = P^{T} S

From where the desired coefficients can be obtained by direct matrix inversion which leads to the classical least-squares matrix inversion method [17],

b = {(P^{T} P)}^{- 1} P^{T} S .

The numerical stability of this solution strongly depends on the chosen polynomial set P. The condition of the Gram matrix G = P^TP = <P,P>, which is the inner product of the polynomials, defines the sensitivity of the solution in Eq. (13) to errors in the data S and is not allowed to be singular. This can be expressed numerically by the spectral condition number

κ (G) = \frac{\max | λ |}{\min | λ |} with λ \in σ (G),

where σ(G) denotes the set of all eigenvalues of G. Condition values near one indicate a well-conditioned matrix while in case of large values, G is said to be ill-conditioned. Numerically, condition numbers are usually estimated, since the exact computation is costly. Determining an estimate for the reciprocal of the condition number 1/κ with values between [0,1] can be performed efficiently using typical solver for systems of linear equations [18]. Here, a value of zero indicates that G is ill-conditioned. However, the direct matrix inversion method was reported to be numerically unstable [5, 19] and needs a tremendous amount of computational effort.

Forsythe [5] pointed out that the solution of Eq. (10) for larger values of polynomial degree M can be greatly simplified by choosing the off-diagonal elements of G_ml(m≠l) to be much smaller than the diagonal elements G_mm. This can be done by selecting a set of polynomials J_m that are orthogonal over the data point set x₁, …, x_N.

Due to the discrete orthogonality condition

G_{m, l} = \sum_{n = 1}^{N} J_{m} (x_{n}, y_{n}) J_{l} (x_{n}, y_{n}) = h_{m} δ_{m, l},

where δ_m,l = 0 for m≠l and δ_m,l = 1 for m = l is the Kronecker delta and

h_{m} = \sum_{n = 1}^{N} {[J_{m} (x_{n}, y_{n})]}^{2}

are the normalization constants, Eq. (10) simplifies to

b_{m} G_{m, m} = c_{m}

where G_m,m = h_m and the coefficients for a best-fit can be found to be simply

b_{m} = \frac{c_{m}}{G_{m, m}} = \frac{c_{m}}{h_{m}} = \frac{\sum_{n = 1}^{N} S_{n} J_{m} (x_{n}, y_{n})}{\sum_{n = 1}^{N} {[J_{m} (x_{n}, y_{n})]}^{2}} .

In case of an orthonormal set of polynomials H_m, with each term normalized by its norm, the associated Gram matrix will have all diagonal elements of unity. Hence, the normalization constants from Eq. (16) will be unity over all orders m and the coefficients can be directly determined by

b_{m} = c_{m} = \sum_{n = 1}^{N} S_{n} H_{m} (x_{n}, y_{n}) .

3. Forbes Q-polynomials

Forbes [7] sag-orthogonal description for general-purpose rotational symmetric aspheres is basically a revised version of Eq. (1) using as well a conic as the base component and an additional set of orthogonal polynomials for high order deformations,

z = f (ρ, ρ_{\max}) = \frac{C ρ^{2}}{1 + \sqrt{1 - C^{2} ρ^{2} (1 + κ)}} + u^{4} \sum_{m = 0}^{M} a_{m}^{c o n} Q_{m}^{c o n} (u^{2}),

including the normalized aperture coordinate u = ρ/ρ_max where ρ_max is the maximum semi-aperture. Forbes described this set to be superior to other choices when it comes to finding a best-fit with a minimum rms sag error to a given surface g(ρ) that can be expressed as

E^{2} (a_{0}^{c o n}, a_{1}^{c o n}, ..., a_{M}^{c o n}) = 〈 {[g (u ρ_{\max}) - \sum_{m = 1}^{M} a_{m}^{c o n} Q_{m}^{c o n} (u^{2})]}^{2} 〉 = minimum,

where the angled brackets represent the weighted average

〈 p (u) 〉 = \frac{\int_{0}^{1} p (u) w (u^{2}) u d u}{\int_{0}^{1} w (u^{2}) u d u},

where w is a weight function that can be chosen freely, e.g. unity for equal weights over the aperture. This will lead, as described in section 2, to the normal equations represented by Eq. (10) with c_l =

〈 g (u ρ_{\max}) u^{4} Q_{m}^{c o n} (u^{2}) 〉

and the Gram matrix

G_{m, l} = 〈 u^{8} Q_{m}^{c o n} (u^{2}) Q_{l}^{c o n} (u^{2}) 〉 = \int_{0}^{1} x^{4} Q_{m}^{c o n} (x) Q_{l}^{c o n} (x) d x,

where x = u². From this, one can define the condition for orthogonality to be

\int_{0}^{1} Q_{m}^{c o n} (x) Q_{l}^{c o n} (x) x^{4} d x = h_{m} δ_{m, l},

in which case the Gram matrix simplifies to G_m,m = h_m with all off-diagonal elements to zero and the normalization constants are given by

\int_{0}^{1} w (x) {[Q_{m}^{c o n} (x)]}^{2} d x = h_{m} .

with the weight function w(x) = x⁴ = u⁸ and the base elements scaled to a maximum value of unity, this condition is fulfilled if the basis is chosen to be a scaled version of the classical Jacobi polynomials [7, 8] P_m⁽^α^,^β⁾(x), which are orthogonal with respect to the weight (1-x) ^α (1 + x) ^β.

Q_{m}^{c o n} (x) = P_{m}^{(0, 4)} (2 x - 1)

Choosing α = 0 and β = 4 complies with the given weight above and yields a standardization where P_m^(0,4)(1) = 1. Since the Jacobi is orthogonal over the interval [-1,1], the additional term in Eq. (26) is needed to adapt its original orthogonality for values from a unit circle, which lie within an interval of [0,1].

The use of recurrence relations allow a recursive evaluation of the basis members to arbitrary high orders based on lower orders in a numerical efficient way. This is especially useful if one wants the mid-spatial frequency components of a surface to be sufficiently described by the surface model. Otherwise, this part will fall into the residuals of the fit where a proper analysis of this component is unfeasible.

The recurrence relation with respect to the degree n is given by [20]

a_{1 n} f_{n + 1} (x) = (a_{2 n} + a_{3 n} x) f_{n} (x) - a_{4 n} f_{n - 1} (x) .

The four corresponding coefficients for the Jacobi polynomials f_n(x) = P_n⁽^α^,^β⁾(x) are well documented in literature. Dividing Eq. (27) by a₁_n results in the standard three term recurrence relation,

f_{n + 1} (x) = (a_{n} + b_{n} x) f_{n} (x) - c_{n} f_{n - 1} (x),

whose new coefficients can be obtained from the former with some calculation. For the special case of α = 0 and β = 4 the coefficients of Eq. (28) are given as [21]

a_{n} = - \frac{(2 n + 5) (n^{2} + 5 n + 10)}{(n + 1) (n + 2) (n + 5)}, b_{n} = \frac{2 (n + 3) (2 n + 5)}{(n + 1) (n + 5)}, c_{n} = \frac{n (n + 4) (n + 3)}{(n + 1) (n + 2) (n + 5)} .

4. Gram-Schmidt process

A set of orthogonal polynomials defined continuously over a certain aperture, loses its orthogonality if used for discrete data sets, especially for lower number of data points N and in case the positions of the data points are irregularly distributed [11, 12]. The integral in the orthogonality condition of Eq. (24) clearly shows that the Q-polynomials are orthogonal for continuous surface sag values over a radial distance between 0 and 1. Therefore, the condition is not fulfilled for a set of discrete surface data points distributed over a circular aperture. To make use of the simplification described in Eq. (18), the polynomial set has to be orthogonalized over the discrete data set using the Gram-Schmidt orthogonalization process which is presented in different ways by a number of classical [13, 14] and modern literature sources [12, 22]. Malacara and DeVore [11] provided a thorough introduction to the orthogonalization process. Mahajan [23] presented a refined version, where the process is extended by an additional normalization step that result in an orthonormal set of polynomials by dividing each orthogonal term by its norm.

The orthogonal polynomials J_m(x, y) and the orthonormal polynomials H_m(x,y) are obtained from the non-orthogonal basis P_m(x, y) in the following manner using abbreviated notation where the arguments of the polynomials are omitted when unnecessary:

J_{1} = P_{1}

H_{m} = \frac{J_{m}}{{‖ J_{m} ‖}_{w}} = \frac{J_{m}}{\sqrt{{〈 J_{m}, J_{m} 〉}_{w}}} = \frac{J_{m}}{{(\sum_{n = 1}^{N} J_{m}^{2} (x_{n}, y_{n}) w (x_{n}, y_{n}))}^{\frac{1}{2}}}

J_{m} = P_{m} - \sum_{k = 1}^{m - 1} D_{m, k} H_{k}

where

D_{m, k} = 〈 P_{m}, H_{k} 〉 = \sum_{n = 1}^{N} w (x_{n}, y_{n}) P_{m} (x_{n}, y_{n}) H_{k} (x_{n}, y_{n})

and w(x_n,y_n) is an arbitrary weight function. Equation (30c) demonstrates the nature of the Gram-Schmidt process, which is a recursive calculation of higher order terms from lower orders by subtracting from each term all components that are parallel to lower order terms. The dependency to the coordinates of the n-th data points (x_n,y_n) in Eqs. (30b) and (30d) indicates that the created polynomial set is only orthonormal for this coordinate set. The newly created orthogonal polynomials differ to the original set. Therefore, fitting to this new set will yield coefficient values a_m that are not connected to the former base set. However, corresponding values for the base coefficients b_m can be retrieved by

D_{m, k} = 〈 P_{m}, H_{k} 〉 = \sum_{n = 1}^{N} w (x_{n}, y_{n}) P_{m} (x_{n}, y_{n}) H_{k} (x_{n}, y_{n})

with a_M = b_M and m = 1, 2, …, M-1 using the conversion matrix

C_{m, l} = \sum_{k = 1}^{K = l - m} D_{l, l - k} C_{l - k, m}

where C_1,1 = 1, m = 2, …, M and l = 1, 2, …, m-1.

5. Numerical performance tests

Various parameters can be investigated to evaluate if the base Q-polynomials are behaving as an orthogonal set and what difference makes the aforementioned orthogonalization or even an orthonormalization. For the evaluation, numerical examples are considered by simulating an aspherical surface using the polynomial expansion of Eq. (1) with the following coefficients

\begin{array}{l} a_{m} = {4 .36653e-7, -2 .21714e-10, -1 .70412e-13, -3 .68093e-17, \\ 8 .94435e-21, 1 .85012e-23, -6 .27043e-27} \end{array}

taken from a typical design of an high-precision aspherical lens. The surface is sampled using an even grid of 101x101 points and cut to a circular aperture leaving a total of N = 7845 data points. The low sampling was chosen to emphasize the effects resulting from the difference between polynomials orthogonal in the discrete and continuous sense.

In an orthogonal set of polynomials the individual terms are not influencing each other and therefore, removing a certain term from the minimization will not affect the coefficients of any other term, making their value independent from the total number of used terms M. Figure 1 demonstrates the change in the coefficients of third and fourth order for different numbers of used polynomial terms. The values are a result of a best-fit in the least-squares sense to the simulated surface using the matrix inversion method from Eq. (13). The remarkable change for different number of terms in case of the base Q^con-polynomials may be regarded as a clear indication that this set, though orthogonal in the continuous case, is not orthogonal for this case of discrete data.

Fig. 1 Selected coefficients of order m = 3 (left) and m = 4 (right) resulting from a best-fit using matrix inversion for different number of polynomials M.

Download Full Size | PDF

At the same time, the coefficients a^ortho of the orthogonal set are completely independent of the number of coefficients, which indicates that the orthogonalization using the procedure from section 3 was successful. The same is true for the orthonormal polynomials though not shown in the graph, since its order of magnitude would scale the plot so that the behavior of the base polynomial cannot clearly be observed.

Another characteristic of orthonormal polynomials is that the mean of each individual term <Q_m> should be zero for all orders m > 0. Therefore, with the first term Q₁ = 1, the total mean of a surface S expanded by these polynomial terms can be estimated by the coefficient a₀. The results in Fig. 2 show that the best estimate can be obtained from the orthonormal set showing the smallest deviations in the order of 10⁻¹⁷. The base polynomials on the other hand show deviations of 15 orders of magnitude larger compared to the results of the orthonormal set. For the orthogonal set, though not as good as the orthonormal set, the deviations are still 13 orders of magnitude smaller compared to the non-orthogonal set resulting in a sufficient estimate.

Fig. 2 Mean value of the individual polynomial terms of order m > 0 for the base set Q^con and its orthogonalized counterparts. Right hand side shows a comparison of the orthogonal terms only.

Download Full Size | PDF

Section 2 presented two different types of solutions to the least-squares problem, the more classical direct matrix inversion method described in Eq. (13) and the special simplified solutions for orthogonal polynomials of Eqs. (18) and (19). The difference between the coefficients resulting from both solutions is plotted in Fig. 3. While the differences for the orthogonal and the orthonormal set are negligible, the differences for the Q^con-polynomials are in the range of 10-100 nm. These values that are within the specification range of modern high precision aspherical surfaces. Therefore, the simple solution is not an appropriate alternative to the classical matrix inversion method for these polynomials.

Fig. 3 Difference between coefficients from simplified solution b_m to coefficients from matrix inversion a_m. Right hand side focuses on the results from the orthogonal terms.

Download Full Size | PDF

However, the matrix inversion method was reported to be numerical unstable. The matrix condition number of the associated Gram matrix and its reciprocal form from chapter 2 can be used as an indicator for the numerical stability of this solution based on the used polynomials. The change of the reciprocal condition number over the number of data points N for different types of polynomials till order m = 16 is plotted in Fig. 4, where values near zero denote an ill-conditioned matrix. Here, the power series expansion of Eq. (1) was included for comparison reasons, demonstrating its obvious ill-condition regardless of the number of data points given. On the other hand, with a constant value of unity, the Gram matrix for the orthonormal polynomials is always well-conditioned independently from N. The base polynomials show a clear dependency where the condition for N < 10³ is similarly ill compared to the power series expansion. For higher values, the condition stays close to ill with a very little improvement of the situation for increasing N. This indicates a certain potential for numerical instability using the matrix inversion method. Surprisingly, the orthogonal polynomials are in a comparable range, indicating the same potential tendency to numerical instability. This illustrates a further advantage of the additional normalization step over the pure orthogonalization.

Fig. 4 Reciprocal condition number of Gram matrix

Download Full Size | PDF

6. Concluding remarks

Fitting measured surface data to model functions for numerical analysis using least-squares minimization is a common task in testing of optical components. Aside from the classical solution of matrix inversion, a simplified solution exists using the inner vector product that is only valid for orthogonal polynomials, as described in chapter 2. This solution is much simpler to realize programmatically, significantly less computationally intensive and therefore, especially suited for embedded systems or real-time applications. However, typical polynomials used in the description of rotationally symmetric aspherical surfaces do not fulfil the orthogonality requirement. As was shown in the previous chapter, even Forbes Q^con-polynomials, an orthogonal set of polynomials over continuous data, loses its orthogonality when applied to discrete data set, especially with smaller numbers of data points. The individual coefficients are influencing each other and show a clear dependency from the number of used polynomial terms. The individual terms are as well not properly normalized and the mean of the surface may not be allowed to be estimated from the coefficient of the first term. Due to the missing orthogonality, the simplified solution may not produce sufficient results. However, a Gram-Schmidt orthonormalization over the discrete data on this basis was used to successfully overcome these drawback yielding an orthonormal equivalent to the basis set for a special given data set. The insignificant difference between the results of both solutions verified the orthonormalized set for the usage of the simple solution.

However, a surface model based on rotationally symmetric descriptions, as the Q-polynomials, is only suitable to characterize surface features with a central axis of symmetry. It is appropriate to model ideal aspherical surfaces and rotational symmetric surface deviations in the mid-spatial frequency range which include common footprints of the rotational machining process of various high precision machining tools. Surface features without a central axis of rotation will not be modelled and will be lost to the residuals of the fit and therefore, cannot be properly characterized. For these cases, orthogonal polynomials like the Jacobi, the Chebyshev of the first kind or the Legendre polynomials, as well as the Zernike polynomials, may be the appropriate choice. The gradient-orthogonal polynomials for freeform surface introduced by Forbes are very promising as well. However, their orthogonality in the slopes makes them less applicable in situations like surface reconstruction by experimental ray tracing that demands a sag-orthogonal basis.

References and links

1. U. Ceyhan, T. Henning, F. Fleischmann, D. Hilbig, and D. Knipp, “Measurements of aberrations of aspherical lenses using experimental ray tracing,” Proc. SPIE 8082, 80821K (2011). [CrossRef]

2. Z. Hosseinimakarem, H. Aryan, A. Davies, and C. Evans, “Considering a Zernike polynomial representation for spatial frequency content of optical surfaces,” in Imaging and Applied Optics, Vol. 1 of 2015, OSA Technical Digest Series (OSA, 2015), paper FT2B.2.

3. E. Goodwin, U. Fuchs, S. Gangadhara, S. Kiontke, V. Smagley, and A. Yates, Design and implementation of a new freeform surface based on Chebyshev Polynomials,” in Imaging and Applied Optics, Vol. 1 of 2015, OSA Technical Digest Series (OSA, 2015), paper FT2B.3.

4. P. Jester, C. Menke, and K. Urban, “B-spline representation of optical surfaces and its accuracy in a ray trace algorithm,” Appl. Opt. 50(6), 822–828 (2011). [CrossRef] [PubMed]

5. G. E. Forsythe, “Generation and use of orthogonal polynomials for data-fitting with digital computers,” J. Soc. Ind. Appl. Math. 5(2), 74–88 (1957). [CrossRef]

6. C. B. Braunecker, R. Hentschel, and H. J. Tiziani, Advanced Optics Using Aspherical Elements (SPIE, 2008), Ch. 2.

7. G. W. Forbes and C. P. Brophy, “Asphere, o asphere, how shall we describe thee?” Proc. SPIE 7100, 710002 (2008). [CrossRef]

8. G. W. Forbes, “Shape specification for axially symmetric optical surfaces,” Opt. Express 15(8), 5218–5226 (2007). [CrossRef] [PubMed]

9. G. W. Forbes, “Fitting freeform shapes with orthogonal bases,” Opt. Express 21(16), 19061–19081 (2013). [CrossRef] [PubMed]

10. D. Ochse, K. Uhlendorf, and L. Reichmann, “Describing freeform surfaces with orthogonal functions,” in Imaging and Applied Optics, Vol. 1 of 2015, OSA Technical Digest Series (OSA, 2015), paper FT2B.4.

11. D. Malacara and S. L. DeVore, “Interferogram Evaluation and Wavefront Fitting,” in Optical Shop Testing, 2nd ed., D. Malacara ed. (Wiley, 1992).

12. V. N. Mahajan, “Zernike Polynomials and Wavefront Fitting,” in Optical Shop Testing, 3rd ed., D. Malacara ed. (Wiley, 2007).

13. G. A. Korn and T. M. Korn, Mathematical Handbook for Scientists and Engineers, (McGraw-Hill, 1968).

14. G. B. Arfken, Mathematical Methods for Physicists, 3rd ed. (Academic, 1985).

15. J. Ye, X. Li, Z. Gao, S. Wang, W. Sun, W. Wang, and Q. Yuan, “Modal wavefront reconstruction over general shaped aperture by numerical orthogonal polynomials,” Opt. Eng. 54(3), 034105 (2015).

16. J. E. Gentle, Numerical Linear Algebra in Statistics (Springer, 1998).

17. J. Y. Wang and D. E. Silva, “Wave-front interpretation with Zernike polynomials,” Appl. Opt. 19(9), 1510–1518 (1980). [CrossRef] [PubMed]

18. E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LAPACK User’s Guide, 3rd ed. (SIAM, 1999).

19. M. P. Rimmer, C. M. King, and D. G. Fox, “Computer program for the analysis of interferometric test data,” Appl. Opt. 11(12), 2790–2796 (1972). [CrossRef] [PubMed]

20. U. W. Hochstrasser, “Orthogonal Polynomials,” in Handbook of Mathematical Functions, M. Abramowitz and I. A. Stegun ed. (Dover, 1978).

21. G. W. Forbes, “Robust and fast computation for the polynomials of optics,” Opt. Express 18(13), 13851–13862 (2010). [CrossRef] [PubMed]

22. R. Upton and B. Ellerbroek, “Gram-Schmidt orthogonalization of the Zernike polynomials on apertures of arbitrary shape,” Opt. Lett. 29(24), 2840–2842 (2004). [CrossRef] [PubMed]

23. V. N. Mahajan and G. M. Dai, “Orthonormal polynomials in wavefront analysis: analytical solution,” J. Opt. Soc. Am. A 24(9), 2994–3016 (2007). [CrossRef] [PubMed]

Fitting discrete aspherical surface sag data using orthonormal polynomials

Abstract

1. Introduction

2. Least-squares minimization with polynomials

3. Forbes Q-polynomials

4. Gram-Schmidt process

5. Numerical performance tests

6. Concluding remarks

References and links

Cited By

Figures (4)

Equations (36)

Optics Express