Robust and precise algorithm for aspheric surfaces characterization by the conic section

A new algorithm for precise characterisation of rotationally symmetric aspheric surfaces by the conic section and polynomial according to the ISO 10110 standard is described. The algorithm uses only the iterative linear least squares. It uses fitting the surface form in a combination with terms containing its spatial derivatives that represent infinitesimal transformations of form. The algorithm reaches sub-nanometre residuals even though the aspheric surface is translated and rotated in the space. he algorithm is computationally robust and an influence of local surface imperfections can be easily reduced by use of a criterion for residuals.


Background
Aspheric surfaces are recently widely used in industry.One of their applications is aspheric lens that often needs its precise characterisation of form.The description of the aspheric surface by the conic section with a polynomial correction is common in ray tracing software and in producer specifications of aspheric lens.The conic section surface fitting with a polynomial correction was addressed by several authors [1][2][3][4][5].Also alternative descriptions of aspheric surfaces were introduced e.g. in [6,7].Nevertheless, a simple and robust algorithm is still needed to evaluate the conic section from measurement data in the ISO 10110-12 form.The coordinate system is shown in the Fig. 1.
The design shape of aspheric surfaces is often described by the z-coordinate as a function of the distance r from z-axis in the form where function c describes the conic section given by function where R is the radius of curvature at the vertex and k is the conic constant (k < -1 hyperbolic, k = -1 parabolic, k > -1 elliptical, k = 0 spherical surfaces).The correction of surface is given by the even-power polynomial with coefficients A.

Radius of curvature
The expansion of (2) shows that the radius of curvature R at the vertex (r is small) can be obtained from linear least squares (i.e.L2-norm) simply using the first coefficient of an even-power polynomial The higher-order terms of Taylor series at the vertex are negligible and the even-power polynomial with the 18 th power is sufficient for reduction of estimation error for R from data with a given range of r.However, high degree of polynomial could introduce numerical errors.They can be reduced by the following way.The linearity of problem allows making the second fit for residuals z ¼ z−z for initial estimate z obtained from q2i as The final estimate of coefficients with reduced numerical error is then thanks to the linearity (additivity) of the problem.The higher power terms are negligible for r close to the vertex.Thus the radius of curvature for rotational paraboloid is obtained as with relatively small error because the first term of expansion of (2) also does not depend on k.The error is below 10 −6 in relative for examples from [3] and the 18 th power polynomial, except the case 1 with relative error 0.002 for R. Thus it is a robust way to evaluate radius of curvature at the vertex.The polynomial with coefficients q 2i also describes the aspheric surface very well for medium precision applications (i.e.1λ flatness of wavefront).Nevertheless, the conic section describes the aspheric form better with less number of coefficients.For example, the Taylor series of (2) for hyperbolic surface with high k converges slowly and thus the even-power polynomial must have more terms for the corresponding precision.

Conic constant
The conic constant k will be obtained by the following way.The initial values are obtained as It corresponds to the parabolic solution from the previous section.The convergence of iterations is worse close to k = -1 because there is a small contribution from the conic section (the terms of expansion for function c).In the next step of algorithm, user selects between the hyperbolic region (k < -1) and the elliptic region (k > -1).If the algorithm output has large errors the second option could be selected automatically.The next values are then R 1 = R 0 and k 1 = -201 or k 1 = -0.5 respectively (The algorithm also works for oblate elliptical surfaces if k 1 is set as a larger positive number.).The value of k 1 for hyperbolic region can be selected closer to the value -1 (e.g.-3 because the most of commercial aspheric lenses have the conic constant above -3).Nevertheless, the value -201 was selected for demonstration purposes.Note that in some cases (e.g. the case 3 from [3]) the maximum radial distance of points r max from axis z is large and it is not possible to calculate the conic section for selected k 1 (e.g.-0.5).In that case the value k 1 -1 is divided by 2 until the k 1 -1 is small enough to calculate the square root in the conic section function c as a real number.Now we calculate two values k 2+ and k 2− (iteration index i = 2) by halving the interval as In the next step, we calculate a pair of differences for all data points using equation and fit them independently by the polynomials p 2j± r 2j using the least squares with selected power larger than r 2 Fig. 1 The sketch of coordinate system used in the ISO 10110 standard and in this paper (e.g. up to r 12 or corresponding to the aspheric lens specifications).The decision between these two fits is based on the lower residual sum of squares.Then the next k and R values are and the iterations are repeated until the residuals are small enough.Then the output coefficients A i are equal to the coefficients p i and the final conic section parameters are k i and R i (Note that the described algorithm also works with negative values of R).The results for five cases from [3] (see also Table 1) are shown in Fig. 2. The convergence is good for all cases and thus the robustness of such algorithm is shown.Nevertheless, the final k values can differ by a few percent from the designed values.
The decision between k i+ and k i-values in the algorithm must be 100% correct.However, an error can occur in some cases.This problem comes from the fact that the combination of conic section with even-power polynomial is underdetermined in parameters within the expected form error. I.e. even thought the relative error of k seems to be large, it cannot be evaluated more precisely from experimental form errors.Nevertheless, the final residuals are in sub-nanometre range even if the initial value k is changed (see Fig. 3).I.e.repeating the procedure with different initial k (in range within the multiple of 2) can improve the results.Also additional conic constants such as e.g.
can be used in each iteration step for the decision of minimal residuals to improve the result (Figs. 2 and 3 show this option).However, it is not necessary because the corresponding residuals are below the uncertainty of measurement that can be carried out.In the case of aspheric lens testing, the known designed value k tbt (that should be calibrated) can be used.Then the initial values can be set e.g.(k 0 + 1) = 0.99(k tbt + 1) and k 1 = k tbt for 1% initial range and the obtained k is then much closer to the lens design value (see Figs. 2 and 3 for case 2).I.e. if the initial parameters are closer to their final values then a less number of iteration steps is needed.

Rotations and translations
The algorithm from previous section and also from e.g.[3] solves the problem where the vertex of aspheric surface is in the origin of coordinates and the surface is not rotated.However, it is not the case for measurement results (3D data of [x, y, z] coordinates) that are generally in an arbitrary coordinate system.This problem can be solved by the following way for relatively small rotations (up to few tens of degrees) and translations (up to few tenths of optical element size).Infinitesimal transformations could be used together (cos α ≈ 1) and thus we can apply linearized substitutions and equation r 2 = x 2 + y 2 into the even-power polynomial equation of r for z.This transformation introduces odd powers of x and y to this polynomial.The initial values of even-power polynomial q 2i are obtained by the linear least squares using equation The additional terms with coefficients g j are used to partially compensate unknown arbitrary transformation in the initial stage.Nevertheless, these coefficients are not used in further calculations.Then the each iteration step consists of the least squares fitting for the following equation where are terms used to represent infinitesimal (derivative) linear and angular transformation contributions to the derived from significant terms solving polynomial expression with substitutions ( 14)-( 16).The even-power polynomial coefficients q2i are taken from the previous iteration step (or they are equal to the initial coefficients q 2i in the beginning).The obtained translation coefficient t z is corrected as t z −α x t x −α y t y (for better convergence of t z ), where coefficients α are total rotation angles from all previous iterations (and are equal to zero at the beginning of algorithm).The resulting translation coefficients t x , t y , t z and rotation angles α x , α y are summed with corresponding translation coefficients t and angles α from the previous iteration.Then the total translation and rotation parameters (i.e. with overline) are used to translate and then to rotate all [x,y,z] values from their original position (to avoid cumulation of numerical errors) and these transformed data points are used in the next iteration step.These iterations are repeated until all transformations do not change with a sufficient precision.Thereafter, the algorithm from previous section is used to calculate aspheric surface parameters such as R and k in correct coordinate system.The rate of convergence of these iterations is superlinear (better for smaller transformations) and faster for higher degree of even-power polynomial.The convergence is slower for larger data sets.Nevertheless, the polynomial degree of 18 (i = 9) is sufficient for iteration convergence of all cases that have been tested (It also includes parameters of various real aspheric lens from different producers.).In the case 1, which is an extreme hyperbolic case with k equal about -171 (not a real lens case), the position error is larger.However, it can be reduced by reduction of the radial range of data or by increasing of polynomial degree.
The good convergence that algorithm reached is shown in Fig. 4. The calculation of position and angles of rotations takes less than one tenth of second on standard computer for data set with a few thousands of points.Thus this algorithm can be used for a real time tracking and displaying of aspheric surface in desired coordinate system.
It should be also noted that in case of large dataset (millions of coordinates), the iterative process could be carried out with some randomly selected part of data (e.g.every thousandth) due to the low roughness of optical surfaces.It speeds up the calculations because the algorithm complexity is given by the linear least squares and thus its calculation time is linearly proportional to the number of data points.The results of such pre-Fig.2 The convergence of relative error of R and k for 5 cases from [3] as a function of the number of iteration index (the polynomial degree was selected the same as the degree of designed A 2i in each individual case).An example for case 2 with the initial value of k within 1% (hollow triangles) Fig. 3 The standard deviation of residuals for 5 cases from [3] and the effect of different initial k on convergence in the case 2 (k = -0.5 is the standard initial value for prolate elliptical surfaces) calculation of transformation and form parameters are used for modification of the initial state of algorithm for subsequent calculation with the full dataset that will describe the form more precisely.

Surface imperfections
The local imperfections in data arise from measurement outliers or effect dust, marks and/or scratches.These imperfections influence the fitting algorithm results although the area of such imperfection is often relatively small.The lens parameters should be evaluated more precisely without an influence of these data (e.g. the focal point of lens).The optical aperture is not so affected if such data from optically not usable areas are removed.The following algorithm with the residual criterion can be used to solve this problem effectively and to keep the robustness of the presented algorithm.
All iteration steps for evaluation of transformation parameters that was described in the previous section will contain the criterion for the correct data.If the residuals of given points will be greater than e.g.10σ (where σ is the standard deviation in z-direction for each iteration step) of least squares residuals then these points will be removed from the data set.
The following example is applied on the case 2 from [3] that is aligned (transformed by the same translations and rotations) as in the previous section.The artificial imperfection is applied to the aspheric form as a difference with the Gaussian form where h is the height of imperfection, s is its width and [x 0 ,y 0 ] are its coordinates.The results for different heights and for s = 0.1 mm, x 0 = 0.5 mm and y 0 = 0.3 mm are shown in Fig. 5.We can clearly see that the error of estimated radius is proportional to the height of such imperfection if it is not filtered out (the filter factor is too large).However, the form is fitted correctly for lower filter factors (such as 10σ in this case).Nevertheless, the filter factor cannot be too small because too many points will be excluded from calculations.Thus some compromise must be made and the corresponding filter factor value is different for different aspheric surfaces and for different Fig. 4 Errors of the radial distances vertex and the total angles of rotation as function of iteration number for 5 cases of surfaces from [3] with initial translations 0.1 mm in all directions X, Y and Z and rotations 0.1 rad and 0.2 rad for axes X and Y respectively Fig. 5 Excluded fraction of points and relative error of estimated radius R as functions of filter factor (sigma multiples) for the case 2 from [3]. (Other parameters such as estimated angle of rotation are not shown because the transition between the correct and incorrect value is similar as for the radius R.) imperfections.In the case of min-max algorithm (e.g. in [8]) the form error is large if 0% of points is excluded.However, it can be significantly smaller if some points are excluded (e.g. about 5%).The parameters obtained from such fitting with a reasonable filter factor can be used to display surface deviations for all points and the standardized form error can be evaluated from these residuals.
Another example of imperfection is measurement noise or roughness of surface.The noise in [x,y,z] data coordinates stops the algorithm convergence.Nevertheless, the resulting error of parameters corresponds to the uncertainty arising from this noise and the algorithm robustness is not affected.
The weighted least squares can be introduced for the conic constant iterative search to deal with the heteroscedasticity of residuals.We can use e.g.weights w corresponding to the projection of orthogonal least squares to the z-direction or to the local curvature of form (i.e.inversely proportional to local radius of curvature) containing the first and the second derivatives of surface analytically calculated from parameters in the previous iteration step.However, the change of weight from w = 1 has a negligible impact on the results.

Results and discussion
The comprehensive test of described algorithm is shown in Fig. 6.An aspheric surface was transformed to a coordinate system unknown for the computer program.Ten artificial imperfections (positive and negative in z-direction) were used and than automatically filtered-out by the fitting algorithm (8.5% of data were removed).The obtained residuals have standard deviation 3.3 pm and maximal error is 11 pm in z-direction (indicated as red colour on the edge of aspheric surface).
The results clearly showed that algorithm allows finding positions of aspheric lenses of unknown form with sub-nanometre accuracy and it was achieved only by least squares fitting.The fitted form for various aspheric surfaces reached residuals also at sub-nanometre level.The non-linear least squares methods (such as in [3]) allows less number of iterations.Nevertheless, linear least squares method, presented here, calculates each step faster and its convergence is robust as it was shown on various examples.All imperfections can be effectively filtered out and precise parameters of lens could be easily obtained.The wide range of algorithm properties such as speed, robustness and effective filtering enables its use in automatic processes.

Conclusions
The algorithm for the evaluation of parameters of rotationally symmetric aspheric surfaces was described and tested.The sub-nanometre precision of fitting was reached with this new robust and fast algorithm.The algorithm also does not need precise estimation of parameters for the initial iteration.It allows ISO 10110 characterization with the conic section parameters of surfaces that are rotated and translated as it is common in the data output from measurement devices and alignment markers on lens are not needed.The additional criterion for residuals can effectively Fig. 6 The fitting errors of algorithm for the case 2 from [3] indicated as colours of points (two views in [mm]).The asphere is translated by 0.1 mm in all directions X, Y and Z and rotated by 0.1 rad and 0.2 rad for axes X and Y respectively.Ten randomly placed Gaussian imperfections were artificially added (with h = 0.05 mm or h = -0.05mm and with s = 0.05 mm).The filter with factor 10σ filtered out these imperfections (holes in the surface) Křen Journal of the European Optical Society-Rapid Publications (2017) 13:11 remove unuseful data to keep results for the main part of aperture unaffected by local imperfections.Moreover, the part of algorithm that finds the correct coordinate system and removes outliers can be used independently on the conic section fitting part.Thus the algorithm will be useful for the optical community for precise characterisation and testing of aspheric lens.