
Asphere is a general term for surfaces that deviate from a sphere^{13}. The phrase ‘aspheric surface’ used herein is specific to a rotationally symmetric aspheric surface. Compared with a spherical surface, which exhibits the same curvature, an aspheric surface exhibits different curvature. Aspheric surfaces have higher degrees of freedom than spherical surfaces, thus allowing them to perform more functions. For example, aspheric surfaces can correct highorder aberrations and improve imaging quality, thus enabling effects that are only possible using multiple spherical mirrors and reducing the size of the optical system. Consequently, optical designers are inclined to use aspheric surfaces in modern optical systems, such as biomedical, lithographic, astronomical optics, and highpower laser systems, when considering system volume and imaging quality^{49}. The design and manufacturing capabilities of aspheric surfaces have gradually improved and highorder aspheric surfaces have become increasingly wellknown owing to the increasing demand for aberration correction in optical systems^{1014}. Sloan et al.^{11} designed a doubleGaussian system and used three eighthorder aspheric surfaces to improve the image quality of an optical system. Yatsu et al.^{12} reduced the zoom lens length of a camera by 30% using only a 10thorder aspheric surface instead of spherical surfaces. Optical systems with aspheric surfaces are continuously being improved owing to the ongoing development of optical design theory^{1526}. The use of aspheric surfaces, particularly highorder surfaces, in optical systems is expected to increase.
The measurement of aspheric surfaces is vital to their manufacture. Although fabrication technology has developed rapidly in recent years^{2735}, advances in fabrication technology have occasionally exceeded measurement capabilities. A wellestablished corollary of the saying, ‘you cannot make it if you cannot measure it’ has been established recently, i.e. ‘If you can measure it, you can indeed make it’. Manufacturing process bottlenecks typically arise from a lack of routine, costeffective, and timely metrology solutions^{36}.
An aspheric surface is primarily measured based on its form and parameters. The surface form is the threedimensional distribution of the surface in the spatial domain. The measurement result of the surface form is a geometric quantity that is typically expressed in terms of the surface height, which is a function of the (x, y) coordinates in units of length. Meanwhile, the aspheric surface parameters (ASPs) are welldefined quantities that can be derived from the surface form. Measurements of the surface form and parameters can be used to assess the aspheric surface quality. According to the international standard ISO 1011012^{37}, the general expression for a rotationally symmetric aspheric surface comprises a conic term and an even term for the power series. The surface form is expressed as shown in Eq. 1.
$$ z\left( h \right) = \frac{{{h^2}}}{{R\left[ {1 + \sqrt {1  (1 + K){{\left( {\dfrac{h}{R}} \right)}^2}} } \right]}} + \sum\limits_{i = 2}^n {{A_{2i}}{h^{2i}}} $$ (1) In Eq. 1, R, K, and A_{2i} are the ASPs; R is the vertex radius of curvature; K is the conic constant; A_{2i} is the highorder aspheric coefficient; h is the surface height perpendicular to the zdirection; and z is the sag of the aspheric surfaces. Fig. 1a shows a schematic illustration of the change in the quadratic aspheric surface with respect to R when K = 0.6. Fig. 1b shows a schematic diagram of the change in the quadratic aspheric surface with respect to K when R = 40. Fig. 1c shows a schematic illustration of the addition of the fourthorder coefficient A_{4} to the aspheric surface based on the quadratic aspheric surface. The length units of the vertical and horizontal axes in Fig. 1 are in millimetres, and the different line types and colours correspond to the various ASPs. Fig. 1 shows that R directly affects the overall curvature of the aspheric surface, the focal length of the aspheric surface, and the curvature near the vertex. Meanwhile, K determines the surface characteristics and type of aspheric surface^{10}. Table 1 shows that K determines whether an aspheric surface is ellipsoid, spherical, paraboloid, or hyperbolic when it is quadratic. K exerts a minimal effect on the curvature near the vertex. The highorder aspheric coefficient A_{2i} is a highorder term added to the quadratic aspheric surface to further adjust the aspheric surface form, which barely affects the curvature near the vertex but significantly affects the curvature at the edges. The aforementioned parameters typically characterise an aspheric surface and are used in the entire design, fabrication, testing, and adjustment of aspherical optical systems.
Fig. 1 Effects of ASPs on curvature: a vertex radius of curvature R, b conic constant K, and c highorder aspheric coefficients A_{2i}.
Range of K Type of aspheric surface $K \in \left( {0, + \infty } \right)$ Short axis ellipsoid $K = 0$ Sphere $K \in \left( {  1,0} \right)$ Long axis ellipsoid $K =  1$ Paraboloid $K \in \left( {  \infty ,  1} \right)$ Hyperboloid Table 1. Effect of conic constant K
Surface form measurements are typically performed in the manufacturing of optical surfaces. Surface form measurement yields the deviation in the spatial domain between the surface under test and the design model, which can effectively guide processing in an optical shop. The peaktovalley (PV) and rootmeansquare (RMS) errors are generally adopted in surface form measurements to assess the quality of a single surface under test and can indicate the processing precision of the surface under test.
However, in an optical shop, assigning the exact final processing precision of each surface using only the PV and RMS errors of an aspherical optical system is challenging. Theoretically, an optical system can meet the requirements as long as the PV and RMS errors of each surface are sufficiently small. However, extremely small PV and RMS errors are rarely required for each surface because such an excessively high processing target will increase the cost significantly. Moreover, aspheric surfaces with excessive PV and RMS errors may form satisfactory systems in practice. The discussion above shows that the surface form measurement of each surface is insufficient for describing the performance of an optical system. Thus, assigning the exact final processing precision target to each surface using only PV and RMS errors is difficult. Hence, another measurement that can reflect the modulation of each light surface should be performed to evaluate the entire system.
In the context of optical system assessments, parameter measurement is an essential step in the manufacture of optical systems as it yields the light modulation deviation (e.g., the focal length change owing to curvature variations, as shown in Fig. 1) between the surface under test and the design model. This allows us to effectively assess whether the surface can function as expected in the optical system. For example, the measurement results of R can provide feedback to the model of an aspherical optical system in commercial optical design software to determine whether the focal length satisfies the requirements before the system is assembled. Moreover, the tolerance analysis of R can offer an exact processing target for aspheric surfaces prior to processing.
The parameter is a concrete form of the surface form in terms of the specific modulation of light. The surface form measurement yields the deviation in the spatial domain and can effectively guide the processing, whereas the parameter measurement can successfully assign the processing target of each surface by assessing the system performance degradation. A routine, costeffective, and timely metrology solution for manufacturing aspheric surfaces is based on a reasonable combination of these approaches.
ASPs are highly demanded in modern advanced optical systems for ensuring good system performance. For example, the R of the primary mirror of the Giant Magellan Telescope is required to be R = 36,000.0 ± 1.0 mm^{38}, the R of the primary mirror of the James Webb Space Telescope is required to be R = 15,879.722 ± 1.0 mm, and the conic constant is required to be K = −0.9967 ± 0.0005 mm^{39}. Therefore, the measurement accuracy should be higher than the abovementioned fabrication accuracy. In this case, the measurement results of the ASPs can offer a vital reference for aspheric processing, thus providing guidance for the following processing steps as well as feedback on possible errors in the processing equipment.
Moreover, Eq. 1 is not the only expression used for the aspheric surfaces. The aspheric surface expressed in Eq. 1 is known as ‘even asphere’ in some commercial optical design software, such as Zemax. An extended version of Eq. 1 with a complete power series is known as ‘odd asphere’, which is expressed in Eq. 2.
$$ z\left( h \right) = \frac{{{h^2}}}{{R\left[ {1 + \sqrt {1  (1 + K){{\left( {\dfrac{h}{R}} \right)}^2}} } \right]}} + \sum\limits_{i = 1}^n {{A_i}{h^i}} $$ (2) The typically performed characterisation of rotationally symmetric aspheric surfaces adds a highorder power series to the conic surface. This representation is simple and provides substantial freedom in optical design. Theoretically, this representation can characterise any symmetric surface with arbitrary accuracy if i is sufficiently large^{40}. Nevertheless, this representation has been reported to be numerically deficient, and the nonorthogonality of the power series introduces considerable difficulty in the leastsquares fitting of ASPs^{40,41}.
Forbes^{40,4244} proposed Q^{con} and Q^{bfs} polynomials based on Jacobi polynomials to address the challenges above. The Q^{con} polynomial was proposed for characterising a strong aspheric surface with large aspheric deviations and represents the aspheric deviations in the zdirection between the strong asphere and its conic base^{45,46}. The Q^{con} surface is generated by combining a conic base and the Q^{con} polynomial with an orthonormal amplitude. The expressions for a Q^{con} surface is shown in Eq. 3.
$$ z\left( h \right) = \frac{{{h^2}}}{{R\left[ {1 + \sqrt {1  (1 + K){{\left( {\dfrac{h}{R}} \right)}^2}} } \right]}} + {w^4}\sum\limits_{i = 1}^n {{B_i}Q_i^{{\text{con}}}\left( {{w^2}} \right)} $$ (3) In Eq. 3, w is the normalised surface height and is defined as w = h/h_{0}, where h_{0} represents the upper limit of h, B_{i} is a coefficient, and $ Q_i^{{\text{con}}} $ is a polynomial term^{37}.
The Q^{bfs} polynomial was proposed for characterising a mild aspheric surface with a constrained slope and represents the aspheric deviations in the normal direction of the bestfit sphere^{45,46}. The bestfit sphere matched the sag of the aspheric surface at the vertices and edges. The Q^{bfs} surface can enhance the manufacturability of an aspheric surface and is generated by combining a spherical base and the Q^{bfs} polynomial, which contains orthonormal derivatives. A Q^{bfs} surface is expressed as shown in Eq. 4.
$$ z\left( h \right) = \frac{{{h^2}}}{{R\left[ {1 + \sqrt {1  {{\left( {\dfrac{h}{R}} \right)}^2}} } \right]}} + \frac{{{w^2}\left[ {1  {w^2}} \right]}}{{\sqrt {1  {{\left( {\dfrac{h}{R}} \right)}^2}} }}\sum\limits_{i = 0}^n {{C_i}Q_i^{{\text{bfs}}}\left( {{w^2}} \right)} $$ (4) In Eq. 4, C_{i} is a coefficient, and $ Q_i^{{\text{bfs}}} $ is a polynomial term^{37}. Moreover, the Q^{bfs} surface can be generated by combining a conic base and the Q^{bfs} polynomial as follows:
$$ z\left( h \right) = \frac{{{h^2}}}{{R\left[ {1 + \sqrt {1  (1 + K){{\left( {\dfrac{h}{R}} \right)}^2}} } \right]}} + \frac{{{w^2}\left[ {1  {w^2}} \right]}}{{\sqrt {1  {{\left( {\dfrac{h}{R}} \right)}^2}} }}\sum\limits_{i = 0}^n {{C_i}Q_i^{{\text{bfs}}}\left( {{w^2}} \right)} $$ (5) The aspheric surface types are listed in Table 2. An aspheric surface is typically generated using a quadratic surface, power series term, or Qtype polynomial. The ASPs include the vertex radius of curvature R, a conic constant K, and higherorder coefficients A_{2i}, A_{i}, B_{i}, and C_{i}. This review primarily focuses on even asphere because it is the most commonly used. In this case, ASPs include vertex radius of curvature R, conic constant K, and highorder coefficient A_{2i}.
Type Basic surface function Additional term Orthogonality Even asphere $ \dfrac{{{h^2}}}{{R\left[ {1 + \sqrt {1  (1 + K){{\left( {\dfrac{h}{R}} \right)}^2}} } \right]}} $ $ \displaystyle\sum\limits_{i = 2}^n {{A_{2i}}{h^{2i}}} $ Nonorthogonal Odd asphere $ \displaystyle\sum\limits_{i = 1}^n {{A_i}{h^i}} $ Nonorthogonal Q^{con} surface $ {w^4}\displaystyle\sum\limits_{i = 0}^n {{B_i}Q_i^{{\text{con}}}\left( {{w^2}} \right)} $ Orthogonal Q^{bfs} surface based on a conic surface $ \frac{{{w^2}\left[ {1  {w^2}} \right]}}{{\sqrt {1  {{\left( {\dfrac{h}{R}} \right)}^2}} }}\displaystyle\sum\limits_{i = 0}^n {{C_i}Q_i^{{\text{bfs}}}\left( {{w^2}} \right)} $ Orthogonal Q^{bfs} surface based on a spherical surface $ \dfrac{{{h^2}}}{{R\left[ {1 + \sqrt {1  {{\left( {\dfrac{h}{R}} \right)}^2}} } \right]}} $ Orthogonal Table 2. Summary of aspheric surface type and representation formulas
Over the past few decades, methods for measuring ASPs have emerged to satisfy the increasing demand for various optical systems. However, the relevant reviews are currently unavailable. This paper reviews various techniques for measuring ASPs, which we hope will contribute to advancements in the fabrication and testing of aspheric optical elements and their practical applications in diverse fields.
The two core issues for a metrology solution are as follows: where data comes from and where data goes. The phrases ‘where data comes from’ and ‘where data goes’ refer to the measurement data source and the data processing approach, respectively. In many research areas, measurement techniques are typically classified based on the measurement data sources, which are directly related to the measuring instrument; therefore, the classification method is highly intuitive and informative. Meanwhile, unlike conventional physical quantities such as length and weight, ASPs cannot be directly quantified by measuring instruments. Therefore, research on ASP measurement is extremely diverse. Classifying ASP measurement methods using measurement data sources is difficult, and no clear classification exists currently.
Therefore, the ASP measurement methods presented herein are classified into two main categories based on the dataprocessing approach: general fitting and centerofcurvaturebased methods. If the measured data are directly used to calculate the ASPs, then the method is classified as the general fitting method. If the measured data are used to position the corresponding centre of curvature, then the method is classified as the centerofcurvaturebased method. In the centerofcurvaturebased method, ASPs are typically calculated based on the distance between the local surface and the corresponding centre of curvature. Subcategories are then classified using a measurement data source.
This paper presents the basic principles, implementation schemes, and exemplary results of the different approaches. A summary and future outlook on techniques for measuring ASPs are provided at the end.

General fitting methods can be classified into the following three types based on the measurement data source: 1) the direct fitting method (surface form), 2) the interferometric method (wavefront aberration), and 3) the geometrical method (surface slope).
The direct fitting method is typically used to measure ASPs and is introduced in Section 2.1. It includes the following three important aspects: absolute surface form measurement, polynomial, and fitting algorithm. It is worth noting that though the surface form measurement is typically measured using an interferometer, the results obtained cannot be applied in the direct fitting method. This is because interferometry is a relative measurement method that can only yield the surface form deviation between the measured and reference surfaces, i.e., it cannot yield the absolute surface form. The technology for measuring the aspheric surface form is relatively mature. Related issues have been discussed by other scholars^{2,4752} and will not be discussed further herein. Section 2.1 introduces the polynomials and fitting algorithms for ASPs.
Accurate measurement of the aspheric surface form is expensive and inefficient. Therefore, other methods have been devised to measure the characteristics of light modulated by a surface instead of directly measuring the aspheric surface form. The phrase ‘characteristics of the light’ refers to a nonnull system wavefront, longitudinal aberration, spot deviation, or other measurable quantity. Complex formula derivations and optical path parameter measurements are typically required prior to ASP fitting. These methods are introduced in Sections 2.2 and 2.3 based on the type of measured characteristics.

In the direct fitting method, a coordinate measuring machine^{53,54}, a swing arm profilometer^{5557}, a Hartmann sensor^{58,59}, a scanning white light interferometer^{60}, or other instruments^{6172} are used to measure the aspheric surface form. Subsequently, the measured data are mathematically fitted to obtain the ASPs. A typically used mathematical fitting process is described as follows^{73,74}:
Determining the best fit for a model function F to a set of measured data F_{n} with N data points typically requires a minimisation process. A model function F can generally be expanded by a polynomial P of degree I as follows:
$$ F = \sum\limits_{i = 0}^I {{P_i}{A_i}} $$ (6) This process can be described as obtaining coefficients A_{i} to minimise the merit function M, which is expressed as
$$ M\left( {{A_0},{A_1} \ldots {A_I}} \right) = \sum\limits_{n = 1}^N {{{\left[ {{F_n}  \sum\limits_{i = 0}^I {{P_i}\left( {{x_n},{y_n}} \right){A_i}} } \right]}^2}} = {\text{minimum}} $$ (7) In Eq. 7, M is a quadratic function of coefficient A_{i}. According to the extreme value theorem, if M is a differentiable function, then the minimum is obtained when the partial derivatives of the function for each coefficient are zero.
$$ \frac{{\partial M}}{{\partial {A_i}}} = 0 $$ (8) This implies that for Eq. 7,
$$ \frac{{\partial M}}{{\partial {A_i}}} = 2\sum\limits_{n = 1}^N {\left[ {{F_n}  \sum\limits_{j = 0}^J {{P_j}({x_n},{y_n}){A_j}} } \right]} (  {P_i}({x_n},{y_n})) = 0 $$ (9) for all i =0,1, 2, …, I, and J = I. Resolving this equation yields
$$ \frac{{\partial M}}{{\partial {A_i}}} = \sum\limits_{j = 0}^J {{A_j}} \sum\limits_{n = 1}^N {{P_i}({x_n},{y_n}){P_j}({x_n},{y_n})}  \sum\limits_{n = 1}^N {{F_n}{P_i}({x_n},{y_n})} = 0 $$ (10) Thus, the following can be obtained at each minimum:
$$ \sum\limits_{j = 0}^J {{A_j}} \sum\limits_{n = 1}^N {{P_i}({x_n},{y_n}){P_j}({x_n},{y_n})} = \sum\limits_{n = 1}^N {{F_n}{P_i}({x_n},{y_n})} $$ (11) Subsequently, the variables below are introduced:
$$ \left\{ {\begin{array}{*{20}{l}} {{c_i} = \displaystyle\sum\limits_{n = 1}^N {{F_n}{P_i}({x_n},{y_n})} } \\ {{G_{i,j}} = \displaystyle\sum\limits_{n = 1}^N {{P_i}({x_n},{y_n}){P_j}({x_n},{y_n})} } \end{array}} \right. $$ (12) Hence, Eq. 11 becomes
$$ \sum\limits_{j = 0}^J {{G_{i,j}}} {A_j} = {c_i} $$ (13) and can be expanded as
$$ \begin{gathered} {G_{00}}{A_0} + {G_{01}}{A_1} + \cdots + {G_{0J}}{A_J} = {c_0} \\ {G_{10}}{A_0} + {G_{11}}{A_1} + \cdots + {G_{1J}}{A_J} = {c_1} \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \\ {G_{J0}}{A_0} + {G_{J1}}{A_1} + \cdots + {G_{JJ}}{A_J} = {c_J} \\ \end{gathered} $$ (14) Eqs. 13, 14 are the normal equations of the leastsquares data fitting problem: If the determinant of G does not vanish, then a unique set of solutions A_{1}, A_{2}, …, A_{i} exists. This equation can be easily obtained using the classical leastsquares matrix inversion method^{75}.
Eq. 6 can be represented in matrix form as
$$\bf F = PA $$ (15) Subsequently, the following can be obtained using the leastsquares matrix inversion method:
$$ {\bf A} = {\left( {{{\bf P}^{\rm T}}{\bf P}} \right)^{  1}}{{\bf P}^{\rm T}}{\bf F} $$ (16) a. Challenges in conventional leastsquare direct fitting method
As shown in Eq. 16, the leastsquares matrix inversion method inverts the Gram matrix G = P^{T}P, which results in numerical instability in the highorder power series fitting process^{76}.
Forsythe^{73} used regression theory to show that the aforementioned leastsquares matrix inversion method is computationally intensive and utilised a typical numerical model to illustrate the strict tolerance of the method under highorder conditions. When N is large and P (x_{n}, y_{n}) is distributed approximately uniformly on (0, 1), one may expect
$$\begin{split} {G_{i,j}}{\text{ = }}&\sum\limits_{n = 1}^N {{P_i}({x_n},{y_n}){P_j}({x_n},{y_n})} \approx N\int_0^1 {{P^i}{P^j}dP }\\=& N\int_0^1 {{P^{i + j}}dP = } \frac{N}{{i + j + 1}}\end{split} $$ (17) In this case, matrix G is N times the matrix [[(i + j + 1) ^{−1}] (i, j = 0, 1, …, k), which is the wellknown principal minor of order k + 1 of the infinite Hilbert matrix.
$$ {\bf H} = \left[ {\begin{array}{*{20}{c}} 1&{\dfrac{1}{2}}&{\dfrac{1}{3}}& \cdots \\ {\dfrac{1}{2}}&{\dfrac{1}{3}}&{\dfrac{1}{4}}& \cdots \\ {\dfrac{1}{3}}&{\dfrac{1}{4}}&{\dfrac{1}{5}}& \cdots \\ \vdots & \vdots & \vdots &{} \end{array}} \right] $$ (18) Solving linear equations involving the minor components of H is difficult. For example, for k = 9, the order of the principal minor H_{10} is 10 and the inverse H_{10}^{−1} contains elements of magnitude 3 × 10^{12}. In this case, a slight error of 10^{−10} in c_{i} results in an error of approximately 300 in A_{i}, indicating that the solution to Eq. 16 is highly sensitive to errors in F_{n} ^{73}. The simplistic fitting process fails when the terms are in excess, owing to the strict tolerance. However, if I is specified within a certain range, the number of terms will be reduced, resulting in the inability of the expression to adequately describe the aspheric surface.
The discussion above indicates that the performance stability of the leastsquares matrix inversion method depends significantly on the Gram matrix G, which is determined by the polynomial set P. The spectral condition number κ can determine the condition of G.
$$ \begin{array}{*{20}{c}} {\kappa ({\bf G}) = \dfrac{{\max \left \mu \right}}{{\min \left \mu \right}}}&{{\text{with}}\; \mu \in \sigma \left( {\bf G} \right)} \end{array}, $$ (19) where σ(G) denotes the set of all eigenvalues of G. Condition values of approximately 1 indicate that G is wellconditioned, whereas high values indicate that G is illconditioned^{74}.
In the fitting process using the leastsquares matrix inversion method and the power series in Table 1, massive coefficients can be obtained when the peak polynomial coefficient increases with the polynomial order. These coefficients are typically several orders of magnitude higher than the sagittal height. This condition may result in a remarkably small number from the subtraction of two large numbers during the sagittal height calculation. Thus, the effective numbers during the calculation decrease because of the limited accuracy of numerical storage in computers. The loss of effective numbers induces numerical instability in the optimisation algorithm used in the calculation. This problem of fitting ASPs can be solved using polynomials and fitting algorithms.
b. Polynomials
An orthogonal polynomial can sufficiently eliminate the numerical instability of the fitting algorithm during the calculation. Forsythe remarked that setting G_{ih} (i ≠ h) to be substantially smaller than the diagonal elements of G_{hh} will significantly simplify the solution to Eq. 14 for large values of polynomial degree I^{73}. A typical practice is to select polynomials that are orthogonal to the dataset as follows: x_{1}, x_{2}, …, x_{N}. Owing to the discrete orthogonality condition, if the polynomial set P_{i} is orthogonal, then
$$ {G_{i,j}} = \sum\limits_{n = 1}^N {{P_i}({x_n},{y_n}){P_j}({x_n},{y_n}) = \left\{ {\begin{array}{*{20}{c}} 0&{i \ne j} \\ {\displaystyle\sum\limits_{n = 1}^N {{{\left[ {{P_j}({x_n},{y_n})} \right]}^2}} }&{i = j} \end{array}} \right.} $$ (20) and the simple form of Eq. 14 can be obtained as follows:
$$ \begin{array}{*{20}{c}} {{G_{00}}{A_0}}&{}&{}&{}& = &{{c_0}} \\ {}&{{G_{11}}{A_1}}&{}&{}& = &{{c_1}} \\ {}&{}&{ \cdots \cdots \cdots }&{}&{}&{} \\ {}&{}&{}&{{G_{JJ}}{A_J}}& = &{{c_J}} \end{array} $$ (21) Thus, the coefficients for the best fit can be easily calculated as follows:
$$ {A_j} = \frac{{{c_j}}}{{{G_{j,j}}}} = \frac{{\displaystyle\sum\limits_{n = 1}^N {{F_n}{P_j}({x_n},{y_n})} }}{{\displaystyle\sum\limits_{n = 1}^N {{{\left[ {{P_j}({x_n},{y_n})} \right]}^2}} }} $$ (22) Instead of a simple additional power series, Forbes^{40,4244} used a nonstandard orthonormal basis in the aspheric expression formula to avoid an illconditioned Gram matrix and proposed three types of function polynomials: Q^{con}, Q^{bfs}, and Q^{bfs} extended. The Q^{con} and Q^{bfs} polynomials are typically used to characterise symmetrical surfaces, as listed in Table 1. Q^{bfs} extended polynomials are used to define complex surface shapes with local protrusions or depressions and can be extended to characterise freeform surfaces. However, Qtype polynomial is defined as orthogonal over continuous data and may not be orthogonal for discrete data. Meanwhile, data obtained from actual tests are typically discrete. Hilbig et al.^{74} proposed the GramSchmidt orthonormalization of these polynomials over a discrete dataset to solve this problem. Wang et al.^{77} expanded the quadratic portion of an aspheric surface into a series of evenorder terms using binomial expansion to generate discrete orthogonal polynomials. The definition of the fitting error is directly related to the inner products; therefore, Cheng et al.^{41} projected polynomials to the vector space based on their relationship with the inner products. Polynomials were analysed using vector analysis methods, where polynomial problems were transformed into vector problems. This method can achieve high accuracy and effectively solve numerical illconditioning problems.
c. Fitting algorithms
The fitting algorithm determines the performance and efficiency of the direct fitting method. Herein, the main list of fitting algorithms applied to ASP fitting is presented. Zhang^{78} (1997) derived aspheric parameter fitting comprehensively and compared various algorithms. The parameter estimation problem was regarded as an optimisation process. A comprehensive discussion of the minimisation criteria and the robustness of the different algorithms was provided. The importance of selecting an appropriate criterion was highlighted as it can affect the accuracy of the estimated parameters, the computation efficiency, and the robustness to predictable or unpredictable errors^{78}. Gugsa^{60} (2005) used the leastsquares method to determine the bestfit conic of aspheric microlenses based on whitelight interferometry, where the Monte Carlo process was used to estimate the measurement uncertainty of the test results. Sun^{79} (2009) utilised the Gaussian–Newton algorithm to fit aspherical curves and surfaces by minimising the vertical distances. This fitting process converged rapidly in the simulated ideal data and data containing random irregularities. In the same year, Chen et al.^{80} discussed the formfitting of rotationally symmetric aspheric surfaces using the leastsquares method and the Levenberg–Marquardt algorithm, which resulted in rapid convergence. ElHayek^{81} (2014) evaluated the following three algorithms based on their performances on simulated datasets: the limited memory Broyden–Fletcher–Goldfarb–Shannon (LBFGS) algorithm, the Levenberg–Marquardt algorithm, and one variant of the iterative closest point algorithm. The LBFGS algorithm showed linear time complexity with respect to the number of data points, executed faster than the Levenberg–Marquardt algorithm, and was significantly faster than the iterative closest point algorithm.
In general, the directfitting method exhibited good versatility. However, the direct use of a nonorthogonal power series for fitting highorder aspheric surfaces causes numerical instability in the calculation. Robustness and efficiency are prioritised in fitting algorithms. A fitting algorithm is susceptible to nonconvergence when the number of fitting points is small. The efficiency of fitting algorithms for addressing largevolume datasets garnered considerable attention. Therefore, researchers have focused on different orthogonal polynomials as well as rapid and accurate fitting algorithms.

Interferometry is an efficient optical method that is widely used to test optical aspheric surfaces^{10, 47, 8288}. This method requires a compensator to reduce the maximum slope of the residual wavefront of the interferometry system; thus, the density of the interference fringes satisfies the Nyquist theorem. The compensator converts the aspheric wavefront into a spherical or plane wavefront, the latter of which interferes with the reference wavefront. Interferograms are analysed to obtain the target information of the aspheric surface under testing (ASUT). Interferometric measurement methods involving a compensator can be classified into null and nonnull interferometry based on the compensator type^{49,8998}.
Null interferometry adopts a null compensator or computergenerated hologram (CGH), which generates a wavefront consistent with the ASUT to compensate for the normal aberration of the aspheric surface, thus allowing a null interferogram to be obtained^{82,83}. This method can achieve high accuracy; however, a null compensator or CGH can only measure a specific aspheric surface, thus significantly increasing the manufacturing costs. The null interferometric measurement system requires the precise positioning of the ASUT at the confocal position of the measurement system, thereby increasing its complexity. In addition, the manufacturing and testing of the null compensator or null CGH is complex, which limits the wide application of null interferometry for ASPs.
Nonnull interferometry is widely used to measure ASPs. As shown in Fig. 2, the typical nonnull interferometry uses a partial compensator (PC) to construct a partial compensation interference system^{49,92,99}. A PC can be used to test a series of aspheric surfaces, thereby reducing the measurement costs and improving the versatility. Unlike null interferometry, nonnull interferometry features a large residual normal aberration in the system wavefront, which is known as the residual wavefront, thus resulting in a retrace error that must be corrected. The residual wavefront is related to the compensation distance, which is defined as the distance between the final surface of the PC and the vertex of the ASUT^{100}. The compensation distance is important and was adopted widely in previous studies.
Yang et al.^{101} proposed a measurement method for R of an aspheric surface using a nonnull interferometry system. As shown in Fig. 3, this method uses the multiconfiguration of a nonnull interferometer for optimisation and can simultaneously yield the R and surface form of the ASUT. An aspheric shift with precisely determined axial displacements was conducted to realise multiple measurements, and a multiconfiguration model was established accordingly. The R and aspheric surface error were entered into the multiconfiguration model as variables, whereas the experimentally tested residual wavefronts were regarded as optimisation objectives. Subsequently, the actual R was retrieved via a simultaneous optimisation process. Theoretically, the accuracy of the results can be improved by increasing the number of measurements. Nevertheless, such improvements typically require additional complex operations and present cumulative errors. The computational efficiency of this method relies on the estimation of the initial value of R and the initial position of the measured aspheric surface. The relative measurement accuracy of this method can exceed 0.025% when measuring a paraboloid mirror with a nominal R of 819.52 mm.
Fig. 3 Schematic illustration of aspheric nonnull interferometry for different axial positions (Reprinted with permission from Ref. 101 Ⓒ Optica Publishing Group)
The measurement of R using the introduced method does not rely on the absolute positioning of the confocal position or the absolute value of the compensation distance, thereby reducing the adjustment difficulty. In fact, this method is applicable to the measurement of R for annular aspheric surfaces. However, its requirements for accurate axial displacement measurements and for fitting rotationally symmetric aberrations in the measured wavefront are high, and neither K nor A_{2i} is measured.
Hao and Hu et al.^{91, 92, 99, 100, 102104} comprehensively developed nonnull interferometry to measure ASPs both theoretically and experimentally.
First, a new concept known as the best compensation distance was defined to assist the measurement^{100,102}. In a nonnull measurement system, the maximum slope of the residual wavefront changes as the aspheric surface propagates along the optical axis. When the maximum slope and density of the interference fringes reaches the minimum value, the compensation distance is defined as the best compensation distance. This position can be identified by observing the density of the interference fringes. The system wavefront obtained at the best compensation distance is the best compensation wavefront.
Subsequently, a mathematical model is established to define the relationship among the ASPs, system wavefront, and best compensation distance. By considering the McLaurin series expansion of Eq. 1, z can be expressed as shown in Eq. 23.
$$\begin{split}& z = \sum\limits_{i = 1}^n {{D_{2i}}{r^{2i}}} {\text{ = }}{D_2} \cdot {r^2} + {D_4} \cdot {r^4} + {D_6} \cdot {r^6} + {D_8} \cdot {r^8} + \sum\limits_{i = 5}^n {{D_{2i}}{r^{2i}}} \\&\quad\quad\quad\quad\quad\quad\quad\quad i = 1,2,3, \cdots ,N \end{split}$$ (23) The relationships between the ASPs and system wavefront are as follows^{103}:
$$ \left\{ \begin{aligned} & \Delta {D_2}{\text{ = }}\frac{1}{{2\left( {R + \Delta R} \right)}}  \frac{1}{{2R}} =  \frac{{\Delta R}}{{2R\left( {R + \Delta R} \right)}} \\& \Delta {D_4}{\text{ = }}\Delta {A_4} + \frac{1}{8} \cdot \frac{{(K + \Delta K + 1)}}{{{{\left( {R + \Delta R} \right)}^3}}}  \frac{1}{8} \cdot \frac{{(K + 1)}}{{{R^3}}} \\& \Delta {D_6}{\text{ = }}\Delta {A_6} + \frac{1}{{16}} \cdot \frac{{{{(K + \Delta K + 1)}^2}}}{{{{\left( {R + \Delta R} \right)}^5}}}  \frac{1}{{16}} \cdot \frac{{{{(K + 1)}^2}}}{{{R^5}}} \\& \Delta {D_8}{\text{ = }}\Delta {A_8} + \frac{5}{{128}} \cdot \frac{{{{(K + \Delta K + 1)}^3}}}{{{{\left( {R + \Delta R} \right)}^7}}}  \frac{5}{{128}} \cdot \frac{{{{(K + 1)}^3}}}{{{R^7}}} \\& \cdots \cdots \\& \Delta {D_{2n}}{\text{ = }}\Delta {A_{2n}} + \frac{{\left( {2n  3} \right)!!}}{{\left( {2n} \right)!!}} \cdot \frac{{{{\left( {K + \Delta K + 1} \right)}^{n  1}}}}{{{{\left( {R + \Delta R} \right)}^{2n  1}}}} \\&\quad\quad\quad \frac{{\left( {2n  3} \right)!!}}{{\left( {2n} \right)!!}} \cdot \frac{{{{\left( {K + 1} \right)}^{n  1}}}}{{{R^{2n  1}}}} \\& \Delta \varphi = \varphi '  \varphi \approx \frac{{4\pi }}{\lambda } \cdot \sum\limits_{i = 1}^n {\Delta {D_{2i}}} \\ \end{aligned} \right. $$ (24) Eq. 24 establishes the relationship between the best compensation wavefront deviation and errors in the ASPs. In this equation set, ΔD_{2i} is the difference between the coefficients D′ and D of the actual and theoretical aspheric surfaces; ΔR, ΔK, and ΔA_{2n} are the errors in ASPs; and Δφ is the deviation between the actual best compensation wavefront φ′ measured in the actual nonnull interferometer and the ideal best compensation wavefront φ. Nevertheless, in this equation set, ΔD_{2i} can be calculated from the last equation while ΔR, ΔK, and ΔA_{2i} (i = 2, 3, and 4) are unknown. A total of i equations and i + 1 unknown quantities are involved. Therefore, the equations are undetermined and have infinite solutions. Similar to the best compensation wavefront, the best compensation distance changes when an error is observed in the ASPs. Another equation can be established by describing the relationship between the deviation of the best compensation distance and the error in the ASPs^{103}. The detailed equation is complex and is omitted herein.
Finally, measurement methods for the best compensation distance remain unknown. Unlike a spherical interferometry system, a nonnull system does not have a cat’s eye or confocal position. Therefore, measuring the best compensation distance, which is key to the entire system, is difficult.
This issue can be mitigated by adding a removable combined aplanat to the PC, as shown in Fig. 4, to locate the cat’s eye position^{100,105}. The layout of this approach is illustrated in Fig. 5. The distance between the PC and the cat’s eye position can be obtained easily based on the aplanat design. After locating the cat’s eye position, the combined aplanat is removed and the aspheric surface propagates along the optical axis on a precision linearity rail to the position where the density of the interference fringes reaches the minimum to calculate the best compensation distance. However, determining the precise location of this position is difficult because the interferogram fringes are similar near the target position. Thus, an iterative optimisation algorithm is utilised to calculate the real best compensation distance and ASPs. The relative measurement accuracies of R and K yielded by this method exceed 0.02% and 2%, respectively, for a conic surface^{104}.
Fig. 4 Layout of removable combined aplanat (Reprinted with permission from Ref. 100 ⓒ Optica Publishing Group)
Fig. 5 Layout of nonnull interferometry with removable combined aplanat (Reprinted from Ref. 104, with the permission of AIP publishing)
However, the versatility of this approach is insufficient. Different aplanats must be designed for various PCs to eliminate spherical aberrations, thus increasing the test cost. The ASUT must traverse a long distance of approximately the R of the ASUT, which is timeconsuming and risky when measuring a large optical element. The accuracy of the guide rail limits the precision of the compensation distance measurement. Theoretically, this approach can be extended to measure A_{2i}; however, further study is required for confirmation.
The laser differential confocal technique was used in another study to measure the best compensation distance^{103}. A schematic illustration of this approach is shown in Fig. 6. The system can be classified into two subsystems: the partialcompensation interferometer system (which is enclosed by the solid line) and the laser differential confocal system (which is enclosed by the dashed line). The real partial compensation interferometer system measures the system wavefront, whereas the laser differential confocal system measures the best compensation distance. The compensation distance is adjusted until the PV value of the residual wavefront reaches the minimum. The wavefront of the system is obtained using a Fizeau interferometer. The PC and ASUT are mounted together to fix the compensation distance. Subsequently, the best compensation distance can be measured after moving the compensator and ASUT to the laser differential confocal system integrally. This approach increases the measurement accuracy of the compensation distance and the versatility of nonnull interferometry, thus allowing highorder aspheric surfaces to be measured, and is suitable for measuring concave and convex aspheric surfaces. The relative accuracies of R, K and A_{2i} can reach 0.025%, 0.095%, and 3.02%, respectively.
Fig. 6 Schematic illustration of nonnull interferometry based on laser differential confocal technique (Reprinted with permission from Ref. 103 ⓒ Optica Publishing Group)
In this method, the PC and ASUT should be fixed and propagated between the interferometer and laser differential confocal system, which is not feasible when testing large optical elements.
In general, the interferometric method can achieve high accuracy in surface form measurements. The information retrieval of this method is based on light wavelength, thus providing superior traceability to the test results. This method requires a compensator to convert the reference wavefront into an aspheric wavefront, and using multiple compensators enable the comprehensive measurement of various convex or concave aspheric surfaces.
The ASPs were calculated by establishing relationships among the ASPs, system wavefront, and compensation distance. Nevertheless, the measurement of the compensation distance does not satisfy the test requirements in various test scenarios. The aplanat approach requires the propagation of the ASUT over long distances. The laser differential confocal technology system is complex, where the ASUT and PC must be fixed and propagated. Implementing these conditions in tests for optical elements with a long R or large diameter is typically difficult. Therefore, a compensation distance measurement method that can satisfy different application scenarios is urgently required.
For interferometry systems, the setup is complicated and the measurement range is relatively small. Furthermore, interferometry cannot correctly retrieve the wavefront when the ASUT presents a height discontinuity greater than onequarter of the light wavelength between two adjacent pixels. Additionally, the 2π phase ambiguity appears in an unwrapped operation^{106}. These conditions limit the applicability of interferometric methods.

Geometrical methods are cost effective and have been widely used to test the surface form of optical elements^{58, 59}, such as via the Hartmann test. Unlike interferometric methods, a geometric optics approach measures the slope from transverse aberrations at some observation planes near the focal plane instead of the optical path difference^{10}, thus avoiding the 2π phase ambiguity problem. Thus, this method has a wide measurement range.

DiazUribe et al. (2000)^{62} proposed a nullscreen method for testing fastaspheric convex surfaces. Fig. 7 shows the mechanism of this method, which uses a cylindrical screen with a set of stripes, dots, drops, or grids. An image is reflected by a perfect surface and a perfect image is created when the dropnull screen is adopted, as shown in Fig. 8a. Each screen is specifically designed to accommodate the ASUT. If the ideal form of the surface to be evaluated is known, then the observed image shown in Fig. 8b indicates the departure from a perfect array if the surface is not ideal. As shown in Fig. 8c, the differences between the theoretical and experimental arrays can be attributed to surface imperfections and can provide relevant information regarding the form of the surface to be analyzed^{71}. As shown in Fig. 9, this method is simpler than other nulltest methods in terms of structure because it does not require any additional optical elements with a specific design to correct the aberrations of the system under test. Slight deviations can be quantified in an image of the grid by measuring the x and ycoordinates at each crossing point. Additionally, alignment can be performed easily because the setup is entirely along the local axes of each segment and the grid or other designed patterns allow the screen to be positioned easily.
Fig. 7 Flatprinted null screen: a fringes and b drops. (Reprinted with permission from Ref. 71 ⓒ Optica Publishing Group)
Fig. 8 Centroids on CCD: a theoretical, b distorted, and c comparison of a and b. (Reprinted with permission from Ref. 71 ⓒ Optica Publishing Group)
Fig. 9 Layout of null screen method (Reprinted from Ref. 10, Copyright ⓒ 2007 by John Wiley & Sons, Inc. All rights reserved)
This method was further developed and can be applied quantitatively to test convex and concave aspheric surface forms^{6472}. Similar to the Hartmann test, this method involves a quantitative evaluation of the nullscreen method via a numerical integration process, such as trapezoidal integration^{10}, and accumulates errors along the integration path. Furthermore, it can be used to qualitatively test large surface departures by visually examining the grid image at the CCD, which is practical during the preliminary manufacturing step. Hence, the ASPs can be fitted when the surface form has been measured.
In addition to the direct fitting approach, the null screen method allows ASPs to be measured directly using novel algorithms. AguirreAguirre et al.(2017)^{71} proposed a randomised algorithm to recover the R, K, and surface forms of their ASUT. In this approach, instead of performing integration and polynomial fitting, a direct and random process was used to determine the surface form. The process involved proposing a new test surface, where the R and K values of the surface are selected randomly, and the null screen was compared with the initially designed one. For the ASUT, R and K, which yielded a null screen that was the most similar to the reference null screen, were specified with the final values. Integration errors were avoided and thus the total error was reduced significantly. The relative accuracies of R and K were 0.4% and 0.5%, respectively.
AguirreAguirre et al. (2018)^{107} extended the method above to calculate a null screen to test fast convex/concave odd aspheric surfaces with highorder coefficients. This approach can extract the target parameters via a simple procedure; hence, a nullscreen test can be conducted without unmounting the optical system from the polishing machine to test the form of the surfaces in the early stages of polishing, which demonstrates the efficiency of this method. The percentage error of this method was shown to be less than 1.3% for the recovery of highorder coefficients A_{i}.

DiazUribe and CornejoRodriguez^{108,109} proposed the geometric method in the 1980s to measure the R and K of quadratic aspheric surfaces. The setup of this method, which comprises a goniometric table and a lens holder with a micrometric screw, is shown in Fig. 10. In this method, the longitudinal aberrations of the normals to the surface and their corresponding angles are obtained on the nodal bench to calculate R and K. The measurement principle is simple and easy to implement. This method is applicable even when the centre of the surface cannot be used but fails when the surface is extremely close to a spherical surface. The accuracy of this method barely satisfies the current test requirements but can be improved using precise auxiliary instruments.
Fig. 10 Layout of nodal bench (Reprinted with permission from Ref. 108 ⓒ Optica Publishing Group)

Wang^{110} proposed a raytracing method to measure the R of aspheric surfaces. The principle of this method is illustrated in Fig. 11. In this method, spotarray rays generated by a laser and grating are projected onto the ASUT, and the positions and directions of the incident and reflected rays are recorded using two shifting receiving planes. A ray tracing method provides the equations for each ray and the coordinates of its intersection point with the surface under testing. Subsequently, measurements can be performed by applying a surfacefitting algorithm to determine the aspheric vertex locations and an optimisation algorithm to calculate the symmetry axis and curvature centre. This method involves a simple structure and can be used to test the R values of most aspheric surfaces. The relative measurement accuracy R exceeds 0.5%.
Fig. 11 Optical system of experimental raytracing method (Reprinted from Ref. 110, Copyright 2004, with permission from Elsevier)
In general, the geometrical method is simple in terms of structure and has a large measurement range. A numerical integration process can be applied to this method to obtain the surface form, and then the ASPs can be fitted. In addition to direct fitting, various algorithms and devices have been proposed to obtain the ASPs directly. A modified method based on a randomised algorithm that avoids integration errors and low spatial resolution problems can yield high performance. The highest relative accuracies of R, K, and A_{i} reached 0.4%, 0.5%, and 1.3%, respectively, when an odd aspheric surface was used as the ASUT.

This review summarises measurement techniques for ASPs, which can be classified into two categories: general fitting and centerofcurvaturebased methods. The general fitting method is the most typically used method for ASP measurements and can be classified into three types: direct fitting, interferometric, and geometric. The capabilities of these techniques are listed in Table 3.
Method/Technique Capabilities Convex Concave R K A_{2i} General fitting method Direct fitting method ◯ ◯ ◯ ◯ ◯ Interferometric method Nonnull interferometry with aplanat × ◯ ◯ ◯ × NonNull interferometry with laser differential confocal technique ◯ ◯ ◯ ◯ ◯ Geometrical method Nodal bench ◯ ◯ ◯ ◯ × Raytracing × ◯ ◯ × × Null screen ◯ ◯ ◯ ◯ ◯ Centerofcurvaturebased method Interferometry × ◯ ◯ ◯ × Laser differential confocal technique × ◯ ◯ × × Table 3. Overview of different measurement techniques for ASPs
The centerofcurvaturebased method is directly compared with the three other submethods of the general fitting method to facilitate discussion. Owing to the gradual maturity of aspheric surface measurement methods, the direct fitting method has been widely used for aspheric surface parameter measurements because of its high accuracy and versatility. The fitting algorithm and polynomials are closely related in the direct fitting method, where the fitting algorithm primarily focuses on efficiency and robustness. The nonorthogonality of the power series in the general expression results in illconditioning problems in the fitting method. Orthogonal polynomials are generated using the Gram–Schmidt process to solve this problem.
Fig. 15 shows that the other three methods offer different advantages in terms of test accuracy, versatility, and cost control. In particular, Fig. 15 only summarises the extreme cases when the accuracy of the corresponding method is the highest and does not represent the performance of all the cases using this method. The interferometric method obtains the ASPs by establishing relationships among the ASPs, system wavefront, and compensation distance. The test results indicate good traceability. Using different compensators, the interferometric method can comprehensively measure various convex or concave aspheric surfaces with high accuracy. Nevertheless, the measurement of the compensation distance using this method does not satisfy the test requirements in various test scenarios. Despite the high accuracy of the interferometric method, the test setup is complicated, and the measurement range is relatively small. The test setup of the geometrical method is simple, and the testing cost is relatively low. The accuracy of the nodal bench and raytracing methods is insufficient and has not been investigated further. The nullscreen method measures ASPs through a surface slope and has a sizeable measurement range. However, data obtained using this method lack traceability. The accuracy of this method depends on the calibration accuracy. The centerofcurvaturebased method can achieve high accuracy when measuring R. However, it features low versatility and cannot be applied for testing highorder coefficients.
Over the recent decades, intensive efforts have been extended to obtain highly accurate measurements. In fact, different efficient approaches have been devised for different ASUTs in various test scenarios. However, the development of advanced optics has presented new challenges. Complex boundary conditions impose strict requirements on optical shop testing. Hence, a test method that comprehensively considers high accuracy, low cost, and versatility should be devised.
Measurement techniques for aspheric surface parameters
 Light: Advanced Manufacturing 4, Article number: (2023)
 Received: 20 September 2022
 Revised: 28 June 2023
 Accepted: 28 June 2023 Published online: 28 July 2023
doi: https://doi.org/10.37188/lam.2023.019
Abstract:
Aspheric surfaces are widely used in advanced optical instruments. Measuring the aspheric surface parameters (ASPs) with high accuracy is vital for manufacturing and aligning optical aspheric surfaces. This paper provides a review of various techniques for measuring ASPs and discusses the advantages/disadvantages of these approaches. The aim of this review is to contribute to advancements in the fabrication and testing of aspheric optical elements and their practical applications in diverse fields.
Research Summary
Parameter measurement: boost of aspheric manufacturing
Aspheric surface is a general term for surfaces that deviate from a spherical surface. Aspheric surfaces are widely used in advanced optical instruments. Measuring the aspheric surface parameters with high accuracy is vital for manufacturing and aligning optical aspheric surfaces. Qun Hao from China’s Beijing institute of technology and colleagues review various techniques for measuring ASPs and discusses the advantages/disadvantages of these approaches. The aim of this review is to contribute to advancements in the fabrication and testing of aspheric optical elements and their practical applications in diverse fields.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article′s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article′s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.