Open Access
Add to BibSonomy
Abstract
Freeform surfaces play an important role in modern optical systems with compactness and better performance. The fabrication tools tend to impart a structured signature on optical surfaces, called ripple errors, during the freeform surface manufacturing process. The description and extraction of ripple errors for freeform surface fabrication and testing have attracted extensive attention. In this paper, we develop a fast and accurate method to describe ripple errors for the large aperture based on Fourier model coupling. The polynomial expression is transformed into Fourier series form and surface errors are reconstructed by frequency feature extraction combining with the least square method. The high accuracy and efficiency of the proposed method for representing and filtering ripple errors consuming little computer memory are demonstrated using real experimental data. The proposed method offers a robust and powerful tool not only suitable for surface error characterization but also for image filtering and analysis.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Freeform surfaces, which are defined as surfaces without rotational symmetric features, have played an important role in illumination [1], imaging [2,3] and high power laser systems [4,5] to improve performance and decrease system weight and size. The higher performance of optical systems put increased complexity and demands on the manufacturing of the optics. The introduction of modern polishing methods, including Magnetorheological Finishing (MRF) [6], Zeeko bonnet polishing [7,8], and Ion Beam Figuring (IBF) [9], replace the traditional grinding for the fabrication of ultra-precision optical elements. However, the existence of periodic ripples called mid-spatial frequency (MSF) errors on the freeform surfaces is inevitable through these deterministic polishing methods. Signatures of these spatially structured errors are tool-specific and causes of MSF errors are the shape of cutting tools, motion control instability, distribution of the motion path, vibration and workpiece deformation, etc. [10,11]. The cut-off frequency for the MSF is defined by the system application and aperture size. The concept of the frequency band of wavefront errors is proposed by Noll [12]. For small apertures with size $D$, the errors of period $D/10$ ∼ $D/40$ are defined as MSF errors [13]. In the study of the National Ignition Facility (NIF), the small scale ripples with the spatial period between 0.12 and 33 mm for the large aperture are defined as MSF errors [14]. Since no interferometer is capable of spanning this entire spatial frequency region, it is decomposed into two spatial regions, PSD1 (33 ∼ 2.5 mm) and PSD2 (2.5 ∼ 0.12 mm) [4], respectively. The small-angle scatter in the transmitted and reflected beam led by MSF errors degrades the image contrast and generates local flare in high-performance optical systems [15,16]. It is necessary to develop a quantitative description of MSF errors for freeform surface fabrication and testing.
It becomes clear that simple peak-valley (PV), root-mean-square (RMS) and roughness measures are no longer adequate to quantify the spatial distribution of MSF errors. The definition of effective surface quality specs thus holds unprecedented challenges. There are different methods of quantifying MSF from measured surface data, such as power spectral density (PSD) [17], area structure-function (ASF) [18], Strehl ratio (SR) [19] and polynomial representations [20]. For PSD analysis, generally, essential to the calculation is low order term removal known as detrending and the use of a window function to minimize spectral leakage of edge discontinuities [21]. The result of characterization depends on the details of the preprocessing steps and is thus not robust to user-choice induced variation. Although extracting the amplitude of frequency is possible through PSD, the information on the anisotropy of MSF texture will be lost. Polynomial representations are to decompose the surface data with the linear combination of polynomial terms out to very high orders. For circular-shaped apertures, Zernike polynomials are a natural basis set with orthogonality. It has been proved that Zernike polynomials can capture different aspects of MSF structures such as type, orientation, period, and their relative phase [22,23]. The generation of very high orders for Zernike polynomials provides a possibility for predicting the degradation of optical performance caused by MSF errors through Zernike fitting coefficients [24]. Q-polynomials is a set of polynomials that are orthogonal in the surface gradient developed by Forbes [25–27]. The fitting process to arbitrary numbers of Q-polynomials terms within circular domains is realized by coupling FFT to a simple explicit form of Gaussian quadrature [28]. This creates promising options for quantifying and filtering MSF structures from measured surface data [20,29]. Radial basis functions (RBF) are using for the representation and tolerance analysis of the surface figure error [30]. It is not adequate to represent MSF structures account of lacking sensitivity to these spatial errors [31,32]. For the large square aperture applied in anamorphic and high power laser optical systems commonly, two-dimensional (2-D) Chebyshev polynomials and 2-D Legendre polynomials constructed from their one-dimensional (1-D) form are orthogonal over the square domain with respect to weight functions [33]. Compared with 1-D Legendre polynomials, 1-D Chebyshev polynomials correspond much better to the familiar Fourier series [34], which demonstrates the possibility of filtering from fitting coefficients. The enormous potential of these two kinds of polynomials in representing MSF errors remains to be explored.
The polynomial representation of MSF errors using the traditional surface fitting method is unavailable on ordinary desktop computers within the limited computational and memory. This paper develops a fast and accurate method to describe ripple errors for the large square aperture consuming little computer memory. Through transforming the polynomial expression into the Fourier series form, surface errors are reconstructed by frequency feature extraction combining with the least square method. Combining with the frequency characteristics of polynomials, the proposed method could be applied to the filtering of wavefront for the designated frequency spectrum. The accuracy and efficiency of the proposed method for representing and filtering ripple errors are demonstrated experimentally.
This paper is organized as follows. Section 2 elaborates on the implementation process of the proposed algorithm based on 2-D Chebyshev polynomials. Section 3 verifies the high accuracy and efficiency of the algorithm for representing MSF errors through several experiments using real experimental data. Section 4 investigates the potential of the proposed method in quantitative filtering. Section 5 discusses and concludes this paper.
2. Algorithm description
After the manufactured surface is measured by interferometry [35], the error surface over the square aperture is obtained by subtracting the ideal surface ${W_{\textrm{ideal}}}(x,\; y)$ without deviation from the real surface ${W_{\textrm{real}}}(x,\; y)$, and is expressed in terms of Cartesian coordinates as
The error surface $W(x,\; y)$ is the input data for all the following descriptions and the domain has been normalized to $x,\; y \in [ - 1,\; 1]$.
For polynomial representations analysis, the error surface is approximated by the sum of a linear combination of basis polynomial functions, which can be expressed as
where ${F_j}(x,\; y)$ represents the $j$th base polynomial, ${a_j}$ is the weighting coefficient, and ${J_{\textrm{total}}}$ is the total number of polynomial terms. The purpose of the surface fitting is to gain weighting coefficients by minimizing the RMS between the original surface and the reconstructed surface. Theoretically, higher fitting accuracy can be realized by increasing the order of polynomials. However, blindly increasing ${J_{\textrm{total}}}$ will confront with computational difficulties, such as consuming enormous computer resources and ill-conditioned solution caused by singular value, through the traditional least square method (LSM). To overcome these obstacles, we propose a fitting algorithm named frequency feature extraction method (FFEM), shown in the following discussion.For the square aperture, there are two alternative polynomials applied directly. Compared with 1-D Legendre polynomials, the intermediate extremes of Chebyshev polynomials are all equal and they correspond much better to the Fourier series [34], which provides the possibility of filtering application. Thus we use 2-D Chebyshev polynomials structured by Chebyshev polynomials in one dimension to represent MSF errors in this work. The recurrence formula of Chebyshev polynomials ${T_n}(x)$ [36] is
2.1 Coordinates transformation
Choosing one of $x$ and $y$ variables for coordinates transformation. For example, we can first denote $x\; = \; \textrm{cos}{{\theta }_x}$, ${{\theta }_x} \in [0,\; \mathrm{\pi }]$. It is easy to deduce Chebyshev polynomials ${T_n}(x)$ after coordinates transformation as
Then the expression of the error surface in Eq. (5) is rewritten as
withNoting that the error surface is expressed as a linear combination of cosine functions and the aim thus turn to obtain the value of weighting functions ${A_n}(y)$ at location $y$.
2.2 Frequency feature extraction
Substituting the Fourier series form of cosine functions into Eq. (7), the expression becomes
Then discretizing ${{\theta }_x}\; = \; j\Delta \; (\; j\; = \; 0,\; 1,\; \ldots ,\; J-1)$, where even $J$ is the number of sampling points in $x$ dimension and the interval $\Delta \; = \; 2\mathrm{\pi }/J$. Through replacing ${{\theta }_x}$ with $j$, ${W_e}({{\theta _x}},\; y)$ is expressed as
To compare with inverse Discrete Fourier Transform (IDFT) directly, Eq. (11) is rewritten as
withIn addition, the condition $J\; \ge \; 2N$ should be satisfied considering the Nyquist sampling theorem. Equation (12) is completely consistent with IDFT in form. Therefore, ${A_n}(y)$ can be calculated by applying the discrete Fourier transform (DFT) to a one-dimensional sampling sequence. Assuming that the one-dimensional sampling sequence ${W_e}(\; j,\; {y_k})$, where location value ${y_k} \in [ - 1,\; 1]$, as known data, one-dimensional spectrum $X(n,\; {y_k})\; (n\; = \; 0,\; 1,\; \ldots ,\; J-1)$ is obtained by applying DFT to ${W_e}(\; j,\; {y_k})$. There is a functional relationship between previous $N$ terms of $X(n,\; {y_k})$ and ${A_n}({y_k})$ as
withParticularly, using FFT instead of DFT can greatly reduce the computing time doubtlessly with the same results.
In addition, sample locations of the error surface measured by interferometers is ordinarily a grid distribution shown in Fig. 1(a). The interval between two adjacent points is equal in two dimensions. After coordinates transformation in our method, equidistant ${{\theta }_x}$ leads to non-uniform $x$ according to the transformation $x\; = \; \textrm{arccos}{{\theta }_x}$, which results in the distribution of sampling points we need is nonuniformly shown in Fig. 1(b). To obtain sampling values corresponding to non-uniform coordinate points, interpolation is exploited here from the raw data. Practice results show that the smooth interpolation method can effectively improve the fitting accuracy. Therefore, we make use of the bicubic interpolation rather than linear interpolation in this work, which has been verified with higher accuracy through experiments.
Fig. 1. Sample locations over the square aperture defined in $x,\; y \in [ - 1,\; 1]$ (equivalent to ${{\theta }_x} \in [0,\; \mathrm{\pi }]$) and $y \in [ - 1,1]$. (a) Uniform grid sampling for the raw data; (b) Non-uniform sampling (sparse in the middle and dense on both sides) used in the FFEM.
In this section, we provide the way to calculate ${A_n}(y)$, which is a combination of base functions shown in Fig. 2. The purpose of solving is transferred to get fitting coefficients.
2.3 Solving coefficients by the least square method
With the estimates derived in Sec. 2.2 for ${A_n}(y)$, the next step is to determine Chebyshev coefficients $a_m^{\; n}$. Based on the least square method, the process is divided into solving $N$ groups of linear equations expressed as
It is worth noting that the selection of ${y_\textrm{1}},\; {y_\textrm{2}},\; \ldots ,\; {y_K}\; $ is not casual. Uniform sampling will cause the solution to failing. One of the methods to overcome this problem adopted in this work shown as
where ${y_k} \in [ - 1,\; 1]$ is the sampling value and ${y_h} \in [0,\; 1]$ belongs to the one-dimensional Halton sequence. This equation includes domain transformation and edge clustering [37], which result in a successful resolution for fitting coefficients.In this section, the objects of the LSM are different from directly using the LSM to the error surface. The right side of linear equations in Eq. (16) is the frequency feature rather than the test data of the error surface. Therefore, unknown variables on the left side only consist of the Chebyshev coefficients corresponding to a single frequency index $n$. By decreasing unknown variables of the LSM through increasing the number of groups of linear equations, it can effectively reduce computing scale and achieve faster computing speed.
2.4 Surface reconstruction
Multiplying all fitting coefficients and the value of 2-D Chebyshev polynomials at sampling points of the error surface and combing them linearly, the fitting error surface is reconstructed directly. However, this reconstruction process is difficult to implement for consuming a great deal of computational time while representing MSF errors out to very high order polynomials. As the last step of the FFEM, a new reconstruction method is proposed here. Assuming that the number of sampling points for the error surface is $Np \times Np$, the coordinates can be represented by column vectors $x\; = \; [{x_\textrm{1}};\; {x_\textrm{2}};\; \ldots ;\; {x_{Np}}]$ and $y\; = \; [{y_1};\; {y_\textrm{2}};\; \ldots ;\; {y_{Np}}]$ because of uniform grid distribution in Fig. 1(a). Then coordinates transformation to $x\; $ is carried out to obtain ${{\boldsymbol \theta }_{\boldsymbol x}}\; = \; [{\theta _\textrm{1}};\; {\theta _\textrm{2}};\; \ldots ;\; {\theta _{Np}}]$. Through the recurrence formula shown in Eq. (3), preceding $M$ terms of Chebyshev polynomials ${T_m}(y)$ at $y$ can be calculated as
According to Eq. (7), the weighting functions ${A_n}(y)$ at coordinates $y\; $ are
where ${\textbf A}$ is the matrix of weighting functions with the size $Np \times N$. Besides, cosine functions $\textrm{cos}n{{\theta }_x}$ at ${{\boldsymbol \theta }_{\boldsymbol x}}$ are arranged asThe reconstructed error surface is
where the matrix $\mathbf{W}$ represents the reconstructed error surface with size $Np \times Np$. The matrix operation coupling frequency feature effectively decrease computational complexity compared with the traditional term by term superposition.There is an explanation that the core idea of the algorithm can be applied to other kinds of polynomials. Another polynomials for the square aperture commonly, 2-D Legendre polynomials, have a more complex calculation process because every term of 1-D Legendre polynomials after coordinates transformation consists of cosine functions with different angular frequencies. Fortunately, by modifying the details of the algorithm FFEM for 2-D Legendre polynomials, it can also be easily carried out. As for Zernike and Q-polynomials for the circular aperture, the process of surface fitting is similar and can refer to the 2-D Chebyshev polynomials applied in the square aperture. We take Zernike polynomials as an example here to illustrate the extension of the algorithm in the circular aperture. The error surface over the circular aperture is generally defined as $W(r,\; {\theta })$ in polar coordinates, where $r \in [0,\; 1]$ and ${\theta } \in [0,\; 2\pi )$. By taking Zernike polynomials to represent the error surface, it is finally described as
The form of this equation is similar to Eq. (9), which means it can be associated with DFT and a similar solution process is adopted. Multi-group one-dimensional sampling sequences corresponding to different values of $r$ are obtained by sampling the error surface at equal angular intervals and fixed radius. By using DFT for each group of sequences, ${A_m}(r)$ and ${B_m}(r)$ correspond to the imaginary part and the real part of the complex frequency spectrum respectively. Then Zernike fitting coefficients are calculated by using the LSM to Eq. (25). Finally, based on fitting coefficients and Eq. (24) and (25), we obtain the reconstructed error surface.
3. Description of ripple errors
In this section, we investigate the accuracy and efficiency of the proposed method for describing ripple errors through several experiments using real data. All experiments are carried out in the same computing environment. The computing platform is a desktop computer with CPU i5-9600K and 32GB of memory, and the mathematical software is MATLAB of version 2020a. To facilitate the following analysis, we illustrate the relationship between the number of 2-D Chebyshev polynomials terms with the maximum spatial frequency of ripple errors. The higher the frequency of spatial errors, the larger the $N$ and $M$ in Eq. (5) needed for fitting, and then the more terms of polynomials are used. Empirically, the more terms are used for fitting, the more high-frequency errors are captured, and then the better the reconstructed error surface is to approximate the original error surface. The premise here is that the number of terms used in fitting does not exceed the number of terms needed to characterize the error surface completely. Otherwise, blindly increasing the number of terms for fitting can not effectively improve the fitting accuracy and even causes a great computational burden. Therefore, it is necessary to estimate the number of terms needed before carrying out the whole representation. The estimation principle is proposed as
3.1 Fitting accuracy
To investigate the fitting accuracy, ripple errors over the large square aperture are represented by the FFEM. Error surfaces in the paper are all tested from mirrors by the 4D dynamic interferometer. The error surface with a $400 \times 400$ mm aperture is plotted in Fig. 3(a), whereas the original data holds $690 \times 690$ pixels. Hence the interval of adjacent pixels is 0.58 mm. According to the resolution of the error surface, the minimum period of ripple errors captured completely by the interferometer is four times of sampling interval (2.32 mm). Therefore, the effective analysis error regions are low-frequency errors and the PSD1 region. To capture high-frequency errors as much as possible, 4 million terms of 2-D Chebyshev polynomials are used for describing this error surface. The minimum error spatial period that can be characterized is $400\mathrm{\pi }/2000 = 1.08$ mm based on Eq. (27). The reconstructed error surface and the fitting residual error surface are presented in Fig. 3(b) and 3(c), respectively, by using the FFEM. All Chebyshev coefficients are shown in Fig. 3(d) while the index $n$ in the $x$ dimension runs left to right (0 to 1999) and the index $m$ in the $y$ dimension increases from top to bottom (0 to 1999) as well. The range of the color bar in that plot is defined artificially to make the value of high-order coefficients display obviously. Noting that the fluctuation of amplitude in Fig. 3(c) has reached the picometer scale. Furthermore, Fig. 3(e) shows the 1-D integral at 0° of the calculated 2-D PSD in the frequency space. The power associated with any given Fourier component of two solid lines does not exceed the dashed line defined by $\textrm{A} \times {f^{\; \; - 1.55}}$, where $\textrm{A}\; = {\; 1.05\; \textrm{nm}^2}\textrm{mm}$ is an empirically derived constant for finished glass optics and $f\; (\textrm{m}{\textrm{m}^{ - 1}}\textrm{)}$ is the spatial frequency [4]. This demonstrates that ripple errors have been effectively suppressed. Figure 3(f) displays an enlarged view of 1-D PSD. According to two nearly coincident curves in Fig. 3(f) and the RMS of the fitting residual error surface at the picometer scale in Fig. 3(c), it can be concluded that our method has capable of describing ripple errors accurately.
Fig. 3. Ripple errors over the large square aperture are represented by 2-D Chebyshev polynomials. (a) The original error surface; (b) The reconstructed error surface based on the FFEM; (c) The fitting residual error surface obtained by subtracting the reconstructed error surface from the original error surface; (d) 2-D Chebyshev coefficients spectrum (The range of the color bar is defined artificially to display obviously); (e) 1-D PSD corresponding to the original error surface, the reconstructed error surface and the NIF standard line; (f) The partial enlargement of 1-D PSD.
To verify the robustness of the algorithm, the error surface containing different frequency components is characterized by 2-D Chebyshev polynomials by using the FFEM. The error surface with a $\textrm{10}.60 \times 10.60$ mm aperture is plotted in Fig. 4(a), whereas the original data holds $480 \times 480$ pixels. Hence the interval of adjacent pixels is 0.022 mm. The minimum period of ripple errors captured completely is 0.088 mm. Therefore, the effective analysis error regions are low-frequency and whole MSF errors. To capture high-frequency errors as much as possible, 2 million terms of 2-D Chebyshev polynomials are used. The minimum error spatial period that can be characterized is $\textrm{10}\textrm{.6}\mathrm{\pi }\textrm{/1414 = 0}\textrm{.023}$ mm. The experiment shown in Fig. 4 is only different from the input error surface data and the number of terms of the experiment in Fig. 3. The RMS values of fitting residual error surface and the enlarged view of 1-D PSD shown in Fig. 4(c) and (f) demonstrate that the fitting accuracy is very high. The proposed method shows the robustness while representing error surfaces with different frequency components.
Fig. 4. Verifying the robustness of the proposed method FFEM. (a) The original error surface; (b) The reconstructed error surface based on the FFEM; (c) The fitting residual error surface obtained by subtracting the reconstructed error surface from the original error surface; (d) 2-D Chebyshev coefficients spectrum (The range of the color bar is defined artificially to display obviously); (e) 1-D PSD corresponding to the original error surface, the reconstructed error surface and the NIF standard line; (f) The partial enlargement of 1-D PSD.
3.2 Fitting efficiency
We investigate the fitting efficiency of the proposed method compared to the traditional least square method in this section. Since the memory of ordinary desktop computers can not afford by directly using the LSM for the full aperture when the number of polynomial terms is more than 100,000, we decompose the full aperture into multiple square sub-apertures for fitting refer to the full aperture is divided into circular sub-apertures for RBF representation [38]. 2-D Chebyshev polynomials are orthogonal to weight functions on each sub-aperture, while the global Chebyshev coefficients can not be calculated. The error surface represented in this experiment holds $720 \times 720$ pixels with an interval of 0.58 mm, which has been through a low-pass filter with a cut-off frequency of $\textrm{1/2}.5\; \textrm{m}{\textrm{m}^{ - 1}}$ from the raw data to remove the interference of high-frequency errors. Thus the fitting residual error surface is only determined by the accuracy of the algorithm rather than high-frequency errors without description. Figure 5(a) illustrates the legend of the sub-aperture least square method (SLSM) on the left and the FFEM on the right. The sampling points of sub-apertures in one-dimension are 144, 120, 90, 80, respectively, and thus the number of sub-apertures is 25, 36, 64, 81. Noting that too small sub-aperture size will seriously reduce the frequency resolution. For comparison purposes, it is necessary to estimate the number of polynomials terms used for the SLSM. Considering the number of polynomial terms for the full aperture is 360,000, the minimum spatial period $T$ is calculated as 2.19 mm based on Eq. (27). According to the sizes of sub-apertures and $T$, it is easy to estimate the number of terms used for the SLSM is 14400, 10000, 5625, 4444, respectively. Empirically, the number of equations is about 1.5 times the number of variables that can achieve higher accuracy in the LSM, and we find the ill-conditioned solution arises when the sampling points of one sub-aperture and the number of polynomial terms are on the same magnitude. By using bicubic interpolation from the metrology data to increase the number of sampling points, this problem has been resolved. The comparison results of two fitting methods are shown in Figs. 5(b)–(d). In Fig. 5(b), the RMS for the SLSM decreases slowly with the reduction of the aperture size, and it is much greater than the RMS for the FFEM. The other two figures show the time and memory consumed by the program of two methods. Both the time consumed and the memory occupied by the FFEM (only consuming 5.26 seconds) are better than the SLSM. It is doubtless that our method gains higher accuracy and higher efficiency at the same time.
Fig. 5. Comparison of accuracy and complexity by two fitting methods. (a) The legend of the SLSM (four gradient colors to represent four aperture sizes) and the FFEM (the red bar represents using 360,000 terms in the fit); (b) The RMS of fitting residual error surfaces based on different fitting methods; (c) Computational time of two methods; (d) The maximum computer memory footprint during the whole process.
4. Quantitative filtering for ripple errors
In this section, we investigate the quantitative filtering for ripple errors based on coefficients partition while Chebyshev coefficients are calculated using the FFEM. Real error surfaces must consist of different forms of MSF structures, which are different in type, orientation, the length scale of the period, the magnitude of different length scales and their relative phase. The raster patterns with only different orientations in Fig. 6(a) and 6(b) are linear sinusoidal patterns across the square aperture and the ring-like pattern in Fig. 6(c) is a radial sinusoidal pattern with the same spatial period. Figures 6(d)–(f) show Chebyshev coefficients from these three MSF structures. The coefficient distribution is determined by types of structures and tends to increase with the enhance of polynomial order in the effective fitting coefficients region. Therefore, we speculate that different error structures can be separated by coefficients partition, and then the purpose of quantitative filtering can be realized.
Fig. 6. Three typical MSF structures are described by 2-D Chebyshev coefficients. (a) Heightmap of vertical sinusoidal raster pattern; (b) Heightmap of tilted sinusoidal raster pattern; (c) Heightmap of radial sinusoidal raster pattern; (d) ∼ (f) Chebyshev coefficients plots from fitting the MSF structures above.
The experiment has been carried out to verify the effectiveness of this filtering method. The error surface tested from the planar element shown in Fig. 7(a) contains low frequency, MSF, and residual high-frequency errors or noise. 2 million terms of 2-D Chebyshev polynomials are used in the fitting process, which means that the minimum period that can be characterized is 0.89 mm. All Chebyshev coefficients in Fig. 7(b) are calculated by the FFEM (the fitting process consumes 98.06 seconds), and the 2-D PSD distribution of the height data is plotted in Fig. 7(c). The restructured error surface in Fig. 7(d) corresponds to Chebyshev coefficients in Fig. 7(e). We remove two coefficients regions in Fig. 7(e), where sector A has a radius of 300 with an angle of $\theta_{1}$ and sector B has a radius of 400 with an angle of $\theta_{2}$, and there are four sectors removed on 2-D PSD in Fig. 7(f). The two sectors with an angle of $2\theta_{2}$ in the vertical direction corresponding to the fitting coefficients of the A region and its radius can be calculated as ${\omega }/(\pi \textrm{D})\; = \; 300/(\pi \times 400)\; \textrm{m}{\textrm{m}^{ - 1}}\; \approx \; 0.24\; \textrm{m}{\textrm{m}^{ - 1}}$. Other two sectors with an angle of $2\theta_{1}$ in the horizontal direction corresponding to the fitting coefficients of the B region and its radius can be calculated as $\textrm{400/(}\mathrm{\pi }\times 400)\; \textrm{m}{\textrm{m}^{ - 1}}\; \approx \; 0.32\; \textrm{m}{\textrm{m}^{ - 1}}$. It is concluded that there is an interesting relationship between Chebyshev coefficients and 2-D PSD. Therefore, we deduce that quantitative filtering can be realized by coefficients partition.
Fig. 7. The corresponding relationship between 2-D Chebyshev coefficients and 2-D PSD. (a) The error surface contains low frequency, MSF and residual high-frequency errors; (b) Chebyshev coefficients from fitting the error surface (The range of the color bar is defined artificially to display obviously); (c) 2-D PSD calculated from the error surface (The range of the color bar is defined artificially as well); (d) The reconstructed error surface (mainly contain inclined ripple errors) from Chebyshev coefficients in (e); (e) Chebyshev coefficients after removing two sectors from (b), where the angle and radius of the A region are $\theta_{1}$ and 300, and the angle and radius of the B region are $\theta_{2}$ and 400; (f) 2-D PSD calculated from the reconstructed error surface and the removed area corresponds to Chebyshev coefficients removed in (e).
We verify the filtering effect based on coefficients partition. The error surface in this experiment is the same as in Fig. 7(a). Chebyshev coefficients in Fig. 7(b) are divided into three parts. The schematic diagram of coefficients partition is shown in Fig. 8(a). As for the PSD1 region, $[{\omega} ]\; \in [37,\; 503]$ means $37\; \le \; \sqrt {{m^\textrm{2}}\textrm{ + }{n^\textrm{2}}} \; \le \; 503$. In Table 1, the number of coefficients for these three parts is about 1000, 200,000 and 1.8 million respectively. The computing time consumed in the process of reconstructing error surfaces is 0.019s, 0.312s and 1.783s, and the memory requirement is about 4M, 8M and 25M. Compared with the time and memory required for solving all fitting coefficients, the time and memory required for reconstructing error surfaces are insignificant. Three reconstructed error surfaces in Figs. 8(b)–(d) are determined by coefficients of three parts. They are low-frequency errors, ripple errors of PSD1, and residual high-frequency errors, respectively. Figures 8(e)–(g) are the 2-D PSD distribution corresponding to these three reconstructed error surfaces. The spatial errors of different frequencies are indeed separated from the PSD view. In Fig. 8(h), we compare the results of three bandpass filtering methods in PSD1 region. The solid blue line is the 1-D PSD of the error surface in Fig. 7(a) as the baseline. The black and green lines represent two filtering methods based on Fourier transform, where the window functions in frequency space are the Hanning window and the error function window, respectively. The band-pass error function window is defined as
Fig. 8. Quantitative filtering by coefficients partition. (a) Chebyshev coefficients from Fig. 7(a) are divided into three parts (low frequency:$\; \sqrt {{m^\textrm{2}}\textrm{ + }{n^\textrm{2}}} \; < \; 37$, PSD1: $\textrm{37 } \le \; \sqrt {{m^\textrm{2}}\textrm{ + }{n^\textrm{2}}} \; \le \; 503$, residual high-frequency: $503\; < \; \sqrt {{m^\textrm{2}}\textrm{ + }{n^\textrm{2}}} $) ; (b) ∼ (d) Reconstructed error surfaces based on coefficients of three coefficient blocks in (a), and they are low frequency, PSD1 and residual high-frequency errors, respectively; (e) ∼ (g) 2-D PSD calculated from three reconstructed surfaces of (b) ∼ (d); (h) The comparison of three bandpass filtering methods.
Table 1. Coefficients number, time and memory required to reconstruct three error surfaces.
5. Discussion and conclusion
Although we have discussed polynomial representations for the PSD1 region of ripple errors, the proposed algorithm FFEM can describe ripple errors of higher frequency as long as the test error surface data with enough resolution. Besides, the core idea of the proposed method is not only suitable for 2-D Chebyshev polynomials but also 2-D Legendre polynomials, Zernike and Q-polynomials. The quantitative filtering application has potential in filtering and image denoising. The proposed method provides a tool to generate polynomials coefficients for ripple errors. These coefficients are possible to evaluate the effect of MSF errors on image quality through the corresponding relationships between the aberrations of the system and the coefficients in the polynomials [34,39]. Besides, the method provides an accurate fitting tool for the optical design method direct mapping [35] to deal with complex freeform surfaces. The goal to set tolerances for fabrication by describing ripple errors under this work will be investigated in future work.
In this paper, we develop a surface fitting method to describe ripple errors for the large square aperture based on Fourier model coupling. The proposed method FFEM has the advantages of high fitting accuracy and high efficiency. The RMS of residual error surfaces reaches the picometer scale and the robustness of the algorithm has been verified while using real experimental data. In addition, quantitative filtering for ripple errors has been proposed. Our method offers a robust and powerful tool not only suitable for surface error characterization but also image filtering and analysis.
Funding
Science Challenge Project (TZ2016006-0502-02).
Disclosures
The authors declare no conflicts of interest.
References
1. R. Wu, L. Yang, Z. Ding, L. Zhao, D. Wang, K. Li, F. Wu, Y. Li, Z. Zheng, and X. Liu, “Precise light control in highly tilted geometry by freeform illumination optics,” Opt. Lett. 44(11), 2887–2890 (2019). [CrossRef]
2. M. Cayrel, “E-ELT optomechanics: overview,” Proc. SPIE 8444, 84441X (2012). [CrossRef]
3. Q. Wang, D. Cheng, Y. Wang, H. Hua, and G. Jin, “Design, tolerance, and fabrication of an optical see-through head-mounted display with free-form surface elements,” Appl. Opt. 52(7), C88–C99 (2013). [CrossRef]
4. J. H. Campbell, R. Hawley-Fedder, C. J. Stolz, J. A. Menapace, M. R. Borden, P. Whitman, J. Yu, M. Runkel, M. Riley, M. Feit, and R. Hackel, “NIF optical materials and fabrication technologies: an overview,” Proc. SPIE 5341, 84–101 (2004). [CrossRef]
5. N. Fleurot, C. Cavailler, and J. L. Bourgade, “The Laser Mégajoule (LMJ) Project dedicated to inertial confinement fusion: Development and construction status,” Fusion Eng. Des. 74(1-4), 147–154 (2005). [CrossRef]
6. J. E. DeGroote, A. E. Marino, J. P. Wilson, A. L. Bishop, J. C. Lambropoulos, and S. D. Jacobs, “Removal rate model for magnetorheological finishing of glass,” Appl. Opt. 46(32), 7927–7941 (2007). [CrossRef]
7. D. D. Walker, A. T. H. Beaucamp, V. Doubrovski, C. Dunn, R. Freeman, G. McCavana, R. Morton, D. Riley, J. Simms, and X. Wei, “New results extending the precessions process to smoothing ground aspheres and producing freeform parts,” Proc. SPIE 5869, 58690E (2005). [CrossRef]
8. C. R. Dunn and D. D. Walker, “Pseudo-random tool paths for CNC sub-aperture polishing and other applications,” Opt. Express 16(23), 18942–18949 (2008). [CrossRef]
9. Y. Lu, X. Xie, L. Zhou, Z. Dai, and G. Chen, “Improve optics fabrication efficiency by using a radio frequency ion beam figuring tool,” Appl. Opt. 56(2), 260–266 (2017). [CrossRef]
10. F. Tinker and K. Xin, “Fabrication of SiC aspheric mirrors with low mid-spatial error,” Proc. SPIE 8837, 88370M (2013). [CrossRef]
11. B. Zhong, H. Z. Huang, X. H. Chen, W. H. Deng, and J. Wang, “Modelling and simulation of mid-spatial-frequency error generation in CCOS,” J. Eur. Opt. Soc.-Rapid Publ. 14(1), 4 (2018). [CrossRef]
12. R. J. Noll, “Effect of mid- and high-spatial frequencies on optical performance,” Opt. Eng. 18(2), 182137 (1979). [CrossRef]
13. H. Gross, ed., Handbook of Optical Systems, Volume 3: Aberration Theory and Correction of Optical Systems (Wiley-VCH, 2007), Chap. Aberrations.
14. W. H. Williams, “NIF large optics metrology software: description and algorithms” Lawrence Livermore National Laboratory Report UCRL-MA-137950, Mar. 16, 2000, p.1–39.
15. R. N. Youngworth and B. D. Stone, “Simple estimates for the effects of mid-spatial-frequency surface errors on image quality,” Appl. Opt. 39(13), 2198–2209 (2000). [CrossRef]
16. J. M. Tamkin and T. D. Milster, “Effects of structured mid-spatial frequency surface errors on image performance,” Appl. Opt. 49(33), 6522–6536 (2010). [CrossRef]
17. T. Pertermann, J. Hartung, M. Beier, M. Trost, S. Schröder, S. Risse, R. Eberhardt, A. Tünnermann, and H. Gross, “Angular resolved power spectral density analysis for improving mirror manufacturing,” Appl. Opt. 29, 8692–8698 (2018). [CrossRef]
18. L. He, C. J. Evans, and A. D. Davies, “Optical surface characterization with the area structure function,” CIRP Ann. 62(1), 539–542 (2013). [CrossRef]
19. M. A. Alonso and G. W. Forbes, “Strehl ratio as the Fourier transform of a probability density of error differences,” Opt. Lett. 41(16), 3735–3738 (2016). [CrossRef]
20. G. W. Forbes, “Never-ending struggles with mid-spatial frequencies,” Proc. SPIE 9525, 95251B (2015). [CrossRef]
21. P. Z. Takacs, “Minimizing Artifacts in Analysis of Surface Statistics,” in Imaging and Applied Optics Technical Papers, OSA Technical Digest (online) (Optical Society of America, 2012), paper OTu1D.5.
22. Z. Hosseinimakarem, A. D. Davies, and C. J. Evans, “Zernike polynomials for mid-spatial frequency representation on optical surfaces,” Proc. SPIE 9961, 99610P (2016). [CrossRef]
23. B. Honarvar Shakibaei and R. Paramesran, “Recursive formula to compute Zernike radial polynomials,” Opt. Lett. 38(14), 2487–2489 (2013). [CrossRef]
24. Z. Hosseinimakarem, A. D. Davies, and C. J. Evans, “Usefulness of orthogonal basis sets for predicting optical performance of wavefronts with mid-spatial frequency error,” Proc. SPIE 10373, 1 (2017). [CrossRef]
25. G. W. Forbes, “Robust, efficient computational methods for axially symmetric optical aspheres,” Opt. Express 18(19), 19700–19712 (2010). [CrossRef]
26. G. W. Forbes, “Characterizing the shape of freeform optics,” Opt. Express 20(3), 2483–2499 (2012). [CrossRef]
27. G. W. Forbes, “Robust and fast computation for the polynomials of optics,” Opt. Express 18(13), 13851–13862 (2010). [CrossRef]
28. G. W. Forbes, “Fitting freeform shapes with orthogonal bases,” Opt. Express 21(16), 19061–19081 (2013). [CrossRef]
29. A. Corumin, “Q-Spectrum and its application for the characterization of mid-spatial frequencies on polished surfaces,” Appl. Opt. 57(19), 5314–5318 (2018). [CrossRef]
30. J. Ni, T. Yang, Y. Liu, D. Cheng, and Y. Wang, “Description and tolerance analysis of freeform surface figure error using specific-probability-distributed Gaussian radial basis functions,” Opt. Express 27(22), 31820–31839 (2019). [CrossRef]
31. J. Stock, A. Broemel, J. Hartung, D. Ochse, and H. Gross, “Description and reimplementation of real freeform surfaces,” Appl. Opt. 56(3), 391–396 (2017). [CrossRef]
32. X. Zhou, X. Huang, J. Bai, L. Zhao, J. Hou, K. Wang, and Y. Hua, “Representation of complex optical surfaces with adaptive radial basis functions,” Appl. Opt. 58(14), 3938–3944 (2019). [CrossRef]
33. J. Ye, Z. Gao, S. Wang, J. Cheng, W. Wang, and W. Sun, “Comparative assessment of orthogonal polynomials for wavefront reconstruction over the square aperture,” J. Opt. Soc. Am. A 31(10), 2304–2311 (2014). [CrossRef]
34. F. Liu, B. M. Robinson, P. J. Reardon, and J. M. Geary, “Analyzing optics test data on rectangular apertures using 2-D Chebyshev polynomials,” Opt. Eng. 50(4), 043609 (2011). [CrossRef]
35. F. Z. Fang, X. D. Zhang, A. Weckenmann, G. X. Zhang, and C. Evans, “Manufacturing and measurement of freeform optics,” CIRP Ann. 62(2), 823–846 (2013). [CrossRef]
36. W. Sun, L. Chen, W. Tuya, Y. He, and R. Zhu, “Analysis of the impacts of horizontal translation and scaling on wavefront approximation coefficients with rectangular pupils for Chebyshev and Legendre polynomials,” J. Opt. Soc. Am. A 30(12), 2539–2546 (2013). [CrossRef]
37. I. Kaya, K. P. Thompson, and J. P. Rolland, “Edge clustered fitting grids for φ-polynomial characterization of freeform optical surfaces,” Opt. Express 19(27), 26962–26974 (2011). [CrossRef]
38. I. Kaya and J. P. Rolland, “Hybrid RBF and local phi-polynomial freeform surfaces,” Adv. Opt. Technol. 2(1), FW2B.4 (2013). [CrossRef]
39. V. N. Mahajan, “Orthonormal aberration polynomials for anamorphic optical imaging systems with rectangular pupils,” Appl. Opt. 49(36), 6924–6929 (2010). [CrossRef]
