1. Introduction

Freeform surfaces offer flexibility in surface definition and aberration compensation in optical system design [1]. Nevertheless, the complexity of surface definition challenges the freeform fabrication process, where small-tool manufacturing techniques are required. During the material removal process, the small tool footprints on the surface create mid-spatial frequency (MSF) errors [2], also known to aspheric manufacturing, which, if not appropriately controlled or mitigated can negatively affect the performance of the optical systems [3–8].

Toward analyzing the impact of MSF surfaces on image performance, one valid approach has its first step to fit measured MSF surfaces to polynomials, such as Zernike polynomials [9–11], 2D-Qs [12,13], and Rapidly Decaying Fourier-like (RDF) polynomials [14,15]. The coefficient maps can then be filtered to reveal unique spatial and angular frequency information for further analysis [8,13,16]. The challenge is that in order to capture the MSF structures, it may be necessary to fit more than a hundred thousand coefficients.

The conventional least-squares (LSQ) fitting method is most commonly employed in this process. LSQ fitting necessitates the computation of the sag table for each individual polynomial term, as well as matrix calculations to identify the coefficients that minimize the sum square of the fitting residuals [17]. This widely used fitting method is not limited to polynomial basis types as long as the sag table for each basis is available.

There are two potential issues brought by the LSQ fitting method when fitting the MSF structures to polynomials. The first challenge is the potential memory overflow caused by solving a large linear algebraic matrix, the size of which is determined by both the number of fitting terms and the quantity of measured data points. One approach to mitigate this challenge is dividing the fitting terms into several smaller groups and performing fittings consecutively. This workaround leads to non-destructive results if the fit is unique, which is often ensured when using orthogonal polynomials such as Zernike polynomials with adequate sampling. The second challenge is the time consumption. Solving linear algebraic matrices can be time-consuming, especially when we fit hundreds of thousands of polynomial terms with a densely measured MSF surface sag table [10].

To address this second most limiting challenge, we have developed a rapid algorithm that leverages both the orthogonality of Zernike polynomials and the Gauss-Legendre quadrature rule for fitting MSF surfaces to Zernike polynomials, bypassing the need for traditional LSQ fitting. The algorithm will be referred to as the GLQ Zernike fitting algorithm. The GLQ Zernike fitting was inspired by a methodology for fitting 2D-Q polynomials without resorting to conventional LSQ fitting [12,13,16]. The differences in the two algorithms relate to the fact that the Zernike polynomials are sag orthogonal, while the 2D-Qs are slope orthogonal. In the broader context of numerical integration techniques, it is noteworthy that alternative approaches like the cubature schemes over the unit disk, as discussed in [18], offer different perspectives, particularly in multi-dimensional integration.

In the remainder of this paper, in Section 2, we will detail the GLQ Zernike fitting algorithm. In Section 3, we will evaluate the algorithm performance in speed and accuracy and compare it to the LSQ fitting before concluding the paper.

2. Fitting method

We shall first rewrite the Zernike polynomials in terms of the Jacobi polynomials, as shown in Section 2.1. The rewrite will allow first to fit the azimuthal components (See Section 2.2), followed by fitting the radial components (See Section 2.3). Central to the GLQ Zernike fitting algorithm is the fact that Jacobi polynomials are orthogonal, which enables the direct calculation of the coefficients. The well-established orthogonality formula will be leveraged together with the Gauss-Legendre quadrature rule to ensure the numerical accuracy of the required integrals. As the quadrature helps to seamlessly transfer integrations to discrete summations, the discrete sampling points will be chosen based on the recurrence matrix of Jacobi. An example of calculating a related matrix is shown in Section 2.4. A discussion of the LSQ fitting method is presented in Section 2.5.

Before delving into the intricacies of the GLQ Zernike fitting algorithm, it is pivotal to understand the foundational representation of freeform optical surfaces. This understanding is rooted in the Zernike polynomials, a key aspect we will explore in Section 2.1.

2.1 Zernike surface

A freeform optical surface can be described as the summation of a set of Zernike polynomials as

2.2 Azimuthal fitting

The azimuthal fitting isolates the $\cos ({m\theta } )$ and $\sin ({m\theta } )$ components and provides the intermediate fitting results with only radial dependencies ${A_m}(u )$ and ${B_m}(u )$ as

Starting with treating the radial position ${u_i}$ as a known parameter while the azimuthal position ${\theta _i}$ is the variable, we evenly sample the azimuthal angle from 0 to 2$\pi $ and set up an array $S_i^{array}$ for each ${u_i}$ as

Then we apply the FFT to $S_i^{array}$, whose real and imaginary parts are exactly the coefficients for $\textrm{the}\; \textrm{cos}({m\theta } )$ and $\sin ({m\theta } )$ bases, respectively.

Please note that we will repeat this FFT operation multiple times with different radial positions ${u_i}$, whose arrangement should follow a specific numerical set, which will be elaborated on in Section 2.3.

2.3 Radial fitting

After getting all the numerical coefficients for each sampled parameter ${u_i}$, the radial component for each m is shown as

The radial fitting leverages the orthogonality of the Jacobi polynomials $P_n^m$, which can be written as

Next, we aim to compute $c_n^m$ from Eq. (14). Continuous integration, a natural choice for this process, is however unfeasible here because the measured data are discrete point clouds, and the fitting process occurs within a digital computing environment. Therefore, we need a solution to handle discrete data effectively while maintaining numerical accuracy. A powerful tool in numerical integration named quadrature is adopted. The general format of the Gauss-Jacobi quadrature rule is given as

Although the Gauss-Jacobi quadrature rule can be leveraged in Zernike fitting once we have the corresponding roots and weights, we choose not to use the Gauss-Jacobi quadrature rule because its root-weight combination is dependent on the factor $\alpha $ and $\beta $, yet $\beta $ varies with the azimuthal order m in our application, which means that for each m order to be fit, a unique set of $({x_1^m,\; x_2^m,\; \ldots ,\; x_K^m;k_1^m,\; k_2^m, \ldots ,k_K^m} )$ needs to be calculated, increasing the computational complexity.

When $\alpha = \beta = 0$, the Jacobi polynomials turns to a particular case, the Legendre polynomials, thus the Gauss-Jacobi quadrature rule becomes the Gauss-Legendre quadrature rule, which is applied in our fitting process. Equation (15) then reduces to

Now the challenge becomes getting the accurate numerical values for all the K roots for the Legendre polynomials ${P_K}(x )$ and the corresponding weights. Here, we bold matrices and put arrows on top of vectors. Fortunately, the root-weight set $[{{x_i},\; {w_i}} ]$ can be calculated with the eigenvalues ${\lambda _i}$ and the corresponding eigenvectors ${\vec{X}_i}$ of ${\boldsymbol J}_{\boldsymbol K}^0$, the tridiagonal recurrence matrix for Jacobi polynomials that satisfy ${\boldsymbol J}_{\boldsymbol K}^0{\vec{X}_i} - {\lambda _i}{\vec{X}_i} = 0$, as [20]

After we replaced the integration according to the Gauss-Legendre quadrature rule, the value of the unknown coefficients in Eq. (14) can be numerically calculated with

Please note that there is a term ${\left( {\frac{{1 + x}}{2}} \right)^{\frac{m}{2}}} = {u^m}$ in Eq. (19), also the calculation for weight ${w_i}$ requires the recurrence of n. The higher m and n order will lead to more calculation complexity, which may come with numerical errors. With that said, the impact of such errors on the overall fitting process is mitigated by the fact that in MSF surface characterization, the higher-order Zernike coefficients typically have significantly smaller amplitudes. Therefore, while the theoretical precision of these high-order terms is affected due to the recurrence relationship, their practical contribution to the shape description of freeform surfaces is often negligible.

The fitting procedure is summarized in Fig. 1. We also published a MATLAB function named “ZernikeLegendreFit” to disseminate the new fitting algorithm [21]. Note that the MATLAB routine calls a prior routine written by Ilhan Kaya that computes the Zernike polynomials recursively [10].

 

Fig. 1. The flowchart of the GLQ Zernike fitting of a freeform surface. The two dashed boxes on the right calculate related roots and weights to use the Gauss-Legendre quadrature rule in the fitting, and these two can be calculated separately even without knowing the measured freeform surface. The measured freeform surface, the highest required fitting terms $N,M$, and the roots ${x_i}$ of the recurrence matrix are used in the azimuthal fitting process, which provides the coefficients of different $\sin (. )$ and $\cos (. )$ bases. Using the Gauss-Legendre quadrature rule, the radial fitting calculates the final coefficients based on Eq. (19).

Download Full Size | PDF

2.4 Tridiagonal recurrence matrix for Jacobi polynomials

As the three-term recursive formula for calculating Jacobi polynomials can be written as

2.5 Least squares fitting

The fundamental principle for LSQ fitting is to solve the least squares problem represented by ${\boldsymbol A}\vec{x} = \vec{b}$. In the context of our application, $\vec{x}$ represents the vector of all the unknown coefficients ${c_i}$, $\vec{b}$ is the vector of measured surface sag $S({\vec{r}} )$, and ${\boldsymbol A}$ is the matrix consisting of calculated values of arbitrary single-indexed Zernike polynomials ${Z_j}({\overrightarrow {{r_i}} } ) $ at the corresponding surface positions $\overrightarrow {{r_i}} $. The matrix equation can be written as

MATLAB also offers iterative-based LSQ functions such as “bicgstab” and “gmres”. Both require the matrix ${\boldsymbol A}$ to be square, meaning the number of polynomials M must equal the number of sample points N on the surface. Properly sampling the surface for MSF fitting requires further exploration and discussion, which is beyond the scope of this paper. In our study, we utilized all data points on the measured surface to ensure robust fitting results.

In conclusion, we find that matrix calculation is the most suitable method for our needs, outperforming other advanced LSQ techniques. Therefore, we will define direct matrix calculation as the “conventional” LSQ fitting method in this research.

3. Evaluation

As introduced in the previous section, the conventional LSQ fitting necessitates the pre-computation of matrix ${\boldsymbol A}$, which includes each individual polynomial value at every surface map position ${Z_j}({{{\vec{r}}_i}} )$. In contrast, the proposed GLQ method still requires pre-computation on each individual polynomial ${P_j}$ value, but the samples are taken only in the radial direction ${u_i}$ rather than from the 2D vector ${\vec{r}_i}$. The ability to adjust the number of radial sample points, K, based on the highest Zernike term, further enhances the efficiency of the proposed fitting algorithm.

3.1 Examples: a synthetic surface and a measured surface

Here, we examine the efficacy of the GLQ Zernike fitting algorithm by applying it to both synthetic and real-measured MSF surfaces. These surfaces are also fitted using the conventional least-squares approach for comparison.

The synthetic MSF surface we chose is a $1024 \times 1024$ matrix defined by $z = 100\sin ({5\pi x} )$ with unitless circular aperture $- 1 < r < 1$. The real-measured MSF surface came as a $925 \times 925$ matrix. Only the pixels within the circular aperture (the inscribed circle of this matrix) of both surfaces have valid sag values. The Zernike coefficients required for fitting are $0 \le m \le M,\; 0 \le n \le N$ and we use the recursive formula to calculate all the Zernike and Jacobi polynomials [22]. It is worth noting that the $m,\; n$ indices we are fitting here are indeed $\tilde{m},\; \tilde{n}$ in Eq. (4), which have been simplified to $m,n $ in Eq. (5).

We consider both the quality and efficiency of each fitting process. The quality is quantified by the root-mean-square (RMS) and peak-to-valley (PV) values of the fitting residual map compared to the original surface. More specifically, the 98% PV (99% highest value – 1% lowest value) and 90% PV are considered to avoid extreme outlier data points [23]. The performed RMS and PV were evaluated within a circular area of normalized radius u ranging from 0 to 0.99. Please note that we do not mask any data or fitting results out when calculating the percentage PV values,

The efficiency, however, is determined by tracking the total time taken for each method. The time consumption, measured using the “tic()” and “toc()” functions in MATLAB, is detailed in these tables under two sections: “Fitting” and “Eval.” The “Fitting” section documents the time spent deriving the coefficients. In contrast, the “Eval” section records the time spent reconstructing the surface based on the coefficients. The evaluation time is not reported for LSQ fitting because the sag maps of each Zernike term are required for the fitting procedure. That is, when LSQ procedure does its coefficient fitting, it calculates the required sag map already, and these sags map can be directly reused in the evaluation process, making the evaluation time negligible.

We have also introduced a time ratio in our analysis. This ratio, which represents the time consumed by the conventional LSQ fitting method divided by the time taken by our proposed GLQ fitting method, has been calculated under two specific conditions: when only calculating the coefficients and when both the coefficients and the fitting error are evaluated.

Both the LSQ and GLQ fittings were done in MATLAB software (2022b) on the same workstation with one AMD Threadripper 3960X processor and 64 GB memory. No parallel computing was allowed during the fitting process. In the LSQ fitting, we divide the whole $({M + 1} )\times ({N + 1} )$ coefficient map into multiple $1 \times 10$ groups and fit them sequentially using matrix inversion; in the GLQ fitting we set $J = 2M + 1,\; K = 2({N + 1} )$ as the azimuthal and radial sampling factors.

In practice, to reduce the impact of aliasing, it may be helpful to first smooth the original image to remove noise before fitting the MSF structure [16]. However, we did not process that here since this evaluation aims to reveal the characteristics of these two algorithms.

The fitting results for the synthetic and the real-measured MSF surfaces are shown in Table 1 and Table 2, respectively, with the residual fitting error maps shown in Fig. 2.

 

Fig. 2. Zernike fitting results. (a) Surface $z = 100\sin ({5\pi x} )$ (left) is fitted to Zernike polynomials with $0 \le m \le 40$; $0 \le n \le 20$ orders. The absolute residual maps from the conventional LSQ Zernike fitting (middle) and the proposed GLQ Zernike fitting (right) are presented with the RMS and PV values listed in their titles. This synthetic surface is unitless. (b) A measured MSF surface named GESPHERE (left) is fitted to Zernike polynomials with the same parameters. This surface’s aperture diameter is 4 mm, and the sag map has its unit in nm.

Download Full Size | PDF

Table 1. Results of fitting a synthetic ${\boldsymbol z} = 100\sin ({5{\boldsymbol \pi x}} )$ MSF surface

View Table | View all tables in this article

Table 2. Results of fitting a real-measured MSF surface

View Table | View all tables in this article

Upon analyzing the fitting residual generated by these two methods for the synthetic surface fitting, shown in Table 1 and Fig. 2(a), we conclude that the proposed GLQ method is superior for this type of surface. This is evidenced by the RMS error from the GLQ fitting being significantly lower than the RMS from the LSQ method. Conversely, for fitting real measured MSF surfaces, as shown in Table 2 and Fig. 2(b), the fitting quality of both methods is comparable, with the RMS and PV values falling within the same order of magnitude. Nevertheless, the speed advantage of GLQ method is consistent across all the fitting experiments. This advantage becomes more pronounced with increasing polynomial order M; and is initially increasing and then decreasing with order N. In fact, specifically for real MSF surfaces, the GLQ method demonstrates a speed improvement of at least several thousand times compared to the conventional LSQ fitting when generating the fitting coefficient maps, which are capable of implying the MSF structures hidden in the surface.

Please note that the two methods (the proposed one and the conventional LSQ fitting) do not usually provide the exact same fitting results or coefficient maps, even if the RMS error of both fittings is comparable. As an example, shown in Fig. 3, we use both methods to fit the sample surface GESPHERE up to $M = 100,\; N = 50$ Zernike terms while the pre- and after-process procedures for these two fittings are kept the same. Although the RMS and PV values from the two methods are comparable, some differences between the two residual maps are noticeable. As shown in Fig. 3(a), the GLQ method fully fits the vertical raster band, whereas the LSQ method still displays a small residual at the location of the band in its residual map. In this specific application, while the LSQ method may be adequate, we conclude that the proposed GLQ method is superior as it provides a cleaner fit.

 

Fig. 3. The results from Zernike fitting the surface GESPHERE with two different methods. $M = 100,\; N = 50.$ (a) The residual maps from fitting surface GESPHERE with the LSQ fitting method (left) and the proposed GLQ fitting (right). Different residual patterns can be observed while the RMS and PV values are comparable. (b) The absolute difference between the two residual maps in (a). (c) The coefficient maps exported from the two methods in (a), respectively, in ${\log _{10}}(. )$ scale. (d) The absolute difference between the two coefficient maps in (c), in ${\log _{10}}(. )$ scale.

Download Full Size | PDF

This comprehensive evaluation substantiates the advantages of the GLQ Zernike fitting algorithm over traditional LSQ fitting. Specifically, we highlight its potential to deliver fitting results with the same accuracy level while significantly reducing the processing time.

Furthermore, another demonstration focusing on the effect of interpolation is discussed in the supplemental material. We also evaluate the aliasing performance of the proposed GLQ fitting algorithm by comparing it with the existing Gauss-Chebyshev fitting method [13], which was originally implemented to fit surfaces to 2D-Q polynomials. Please refer to the supplemental material (Supplement 1) for more information.

3.2 Generalization beyond the Zernike polynomial

This proposed method employs the Gauss-Legendre quadrature rule to convert integration into discrete summation, a process which is succinctly described Eq. (16). Due to its simplicity, this fitting algorithm can be generalized to other types of polynomials as long as the orthogonality relationship of those polynomials are maintained.

For example, to fit a signal to Chebyshev polynomials shown as $f(x )= \mathop \sum \limits_n {t_n}{T_n}(x )$ on $[{ - 1,\; 1} ]$, as the orthogonality relationship of Chebyshev polynomial is shown as

Combining Eq. (29) and (30), the Chebyshev coefficients can be calculated as

When applying the proposed GLQ fitting method to RDF $({{N_b} = 1} )$ polynomials [14,15], which have a well-defined radial component with clear orthogonality as indicated in Eq. (30) and Eq. (27a) in [14], respectively, we encounter a specific challenge. Although we can fit MSF surfaces to RDF using GLQ, the lack of a rapid and accurate calculation method for $\mathrm{\Re }_n^m(r )$ diminishes the speed advantage reported in fitting Zernike polynomials.

4. Conclusion

In a manner akin to the approach in the previous paper [13], which implemented a robust method for fitting MSF surfaces to 2D-Q polynomials, our study is motivated to fit Zernike polynomials, which are more commonly used in various contexts. The fitting algorithm for Zernike polynomials developed in this work, which leverages the Gauss-Legendre quadrature rule, has proven to maintain equivalent fitting accuracy while significantly increasing fitting speed. It is also noteworthy that [13] discusses an alternative approach for fitting Zernike polynomials, which are compared with the proposed Gauss-Legendre fitting in the supplemental material with respect to fitting aliasing, contributing to the diversity of methods available for surface characterization.

Our study demonstrated a quadrature fitting algorithm capable of computing Zernike coefficients without the necessity of evaluating each individual polynomial term at every map position. This process increases the speed by over three orders of magnitude (i.e., more than 1000 times) compared to the conventional least-squares method, when only the coefficient map is needed without a subsequent surface reconstruction. Even when surface reconstruction is required – a process involving the calculation of the surface sag map based on the fitted coefficients – the Zernike GLQ fitting algorithm remains significantly more efficient, approximately three times faster than the conventional approach, highlighting the method’s robust efficiency in diverse application scenarios.

A point to note is that the proposed GLQ fitting algorithm does have limitations, particularly its dependence on interpolation. Since interpolation aims to generate a smooth surface or curve from surrounding data points, it is important to note that this method is unsuitable for signals with large areas or sections removed within their aperture, as they may not be accurately reconstructed.