1. Introduction

Freeform optics play established roles for progressive ophthalmic lenses, in conformal and non-imaging applications, and across a diversity of display and projection systems. Several conventions have been developed and implemented for describing their surface shapes in commercial optical design codes. Perhaps the most standard option is a simple polynomial added to a conic section, but a range of others exists, see e.g [1–4]. While diamond turning machines can deliver parts with the desired figure and finish in one step for some applications (e.g. infrared and illumination), higher precision contexts generally require grind-and-polish processes. For this latter category of molds and parts, a specially tailored gradient-orthogonal basis was introduced to characterize shape [5]. While this basis was constructed to be robust and efficient and to facilitate the assessment of manufacturability, it was recently shown to deliver another benefit. In particular, the associated coefficients were found to be especially effective degrees of freedom during design [6]. Although they span the same space of surface shapes as other common representations, solutions with superior optical performance were discovered when optimizing in terms of this new description. Together with its efficiency, the gradient-orthogonal basis therefore offers significant promise for applications beyond those driven by consideration of “grind-and-polish manufacturability”.

Effective means for interconverting alternative descriptions of shape are valuable because no single convention can meet all needs. One of the advantages of an orthogonal basis is that fitting any particular shape in such a basis can often be performed more efficiently by exploiting orthogonality. Since the values of the associated fit coefficients are independent of the particular subset of the basis elements that is in use and follow directly from simple inner products, the resulting fitting process is sometimes referred to as an “orthogonal projection” onto the basis. A projection of this sort is reported here for the gradient-orthogonal basis of [5]. The process is constructed to avoid the need for determining derivatives and is shown to open the possibility of new applications of these ideas for the analysis of mid-spatial frequency structure (MSF) on manufactured parts.

Modern sub-aperture tools lack the automatic smoothing of traditional spherical fabrication processes. This has given unprecedented importance to both quantifying and tolerancing MSF, but these steps continue to be problematic. Although it may seem intuitive to use Fourier methods to determine a power spectral density (PSD), the lack of data beyond the part’s aperture causes the well-known phenomenon of “spectral leakage”. Various window functions have been introduced to control this corruption, but these inherently mask edge effects on the part that are oftentimes of critical interest. “Detrending” is another standard step for reducing the impact of spectral leakage. This involves subtracting a low-order polynomial from the data to help to suppress the artificial discontinuity at the aperture’s edge before applying a window and Fourier transform. The arbitrary aspects of these steps mean that interpretations of the end results are troubled with questions. It is the finite aperture that generates these complications so it is natural to seek a Fourier-like decomposition that is matched to the aperture at its foundation. The robustness and efficiency of the general sorts of processes described below create such an option for fitting manufactured part shapes entirely with polynomials. The coefficients in the resulting fit are shown to provide a new spectrum that enables familiar types of spatial frequency filtering operations and a promising quantification of MSF.

In keeping with [5], the sag of the freeform surfaces considered in this work is expressed in terms of cylindrical polar coordinates as $z=f\left(\rho ,\theta \right)$ where, with $u=\rho /{\rho }_{\max }$,

When fitting the parameters in Eq. (1.1) to some other representation of sag, say $z=S\left(\rho ,\theta \right)$ where the function $S$ is regarded as a black box, the objective is to determine values for $a_n^m$, $b_n^m$, and $c$ to match $f\left(\rho ,\theta \right)$ to $S\left(\rho ,\theta \right)$. The first step is to ensure the best-fit sphere has the required sag values at $\rho =0$ and at $\rho ={\rho }_{\max }$. Because $f(0,\theta )$ evaluates to zero, the alternative sag description is given a constant bias of $S_0:=-S\left(0,\theta \right)$, where “$:=$” means “is defined to be”. This offset in $z$ simply translates the surface without changing its shape. If averaging over angle is denoted by

Fitting to a linear combination of, say $E$, basis elements generally requires the solution of a set of linear equations involving an $E\times E$ Gram matrix [7]. For an orthogonal basis, however, this process simplifies significantly. To exploit orthogonality for Eq. (1.7), the fit must be performed over the entire enclosing disc (i.e. $0<u<1$ and $0<\theta <2\pi$) even if the surface’s clear aperture is only a subset of that disc. The simplification starts, for this case, by splitting the process into two phases. In particular, rather than solving for all of the fitting coefficients in a single step, it is significantly more efficient to break the operation into separate stages of azimuthal and radial conversions. First, for any specific value of $u$, it is straightforward to determine $A_m(u)$ and $B_m(u)$ so that

The Fourier series methods used for Eq. (1.8) are summarized in Sec. 2. That brief review includes a discussion of the impact of aliasing and, importantly, it gives an idea of an analogous issue in working with Eqs. (1.9) and (1.10) to complete the projection. A powerful method for solving Eq. (1.9) for $a_n^0$ has already been developed [8] and its end result is summarized below in Sec.3. As discussed in Sec. 2, standard Fourier methods involve a particularly simple variant of Gaussian quadrature. Similarly, as pointed out in Sec. 4, the optimal solutions for Eqs. (1.9) and (1.10) involve what is known as Chebyshev-Gauss quadrature [9]. Demonstrations are presented in Secs. 5 and 6 before the Concluding Remarks of Sec. 7. Readers who seek further inspiration to take on the mathematical details may wish to jump first to the figures of Sec. 6.

2. Fourier methods for the azimuthal fit

Although the Fourier series is elementary, some of its aspects offer helpful insights for the sections that follow. For clearer access to the points of interest, I initially express the results in terms of complex exponentials; they are later recast in terms of explicitly real-valued functions. The Fourier series for a function of period $2\pi$ can be written as

When the coefficients, i.e. ${\ell }_k$, are chosen to minimize the mean squared modulus of the difference between the left- and right-hand sides of Eq. (2.1), it follows from Eq. (2.2) that the associated Gram matrix is an identity. These coefficients are therefore given directly by

Notice that the dominant computation for the fit of Eq. (2.1) by using Eq. (2.3) lies in numerically approximating that integral for all the required values of $k$. Valuable results emerge if a uniformly distributed and weighted sum is used for that approximation. For example, consider the case of sampling at $\theta =j\Delta$ for $j=1,2,3,\mathrm{...}J$ where $\Delta =2\pi /J$ (with, say, $\alpha =\Delta /2$ in mind). The integral in Eq. (2.3) is then approximated by

The mixing of terms in Eq. (2.6) is referred to as aliasing. For band-limited functions, i.e. where ${\ell }_k$ vanishes whenever $\left|k\right|>B$ for some $B$, it follows that ${\tilde{\ell }}_k$ is precisely equal to ${\ell }_k$ over all of the non-zero band provided $J>2B$. That is, if we sample sufficiently densely, all the terms but $n=0$ vanish in Eq. (2.6) for $\left|k\right|\le B$. This requirement that $J>2B$ is the familiar Nyquist sampling condition. In most cases of interest to us, $\left|{\ell }_k\right|$ decays rapidly (oftentimes exponentially) as $\left|k\right|$ grows, so the functions are effectively band limited and aliasing has negligible impact when sufficiently dense sampling is used, i.e. $J$ is sufficiently large. The second objective in what follows is to find analogous summation-based processes that reduce the computational burden of the steps involved in the final stage of projecting to$Q_n^m$, i.e. for working with Eqs. (1.9) and (1.10).

Before proceeding, it is helpful to notice that when $L(\theta )$ is real valued, it follows from Eq. (2.3) that ${\ell }_{-k}={\ell }_k{}^{\ast }$. Now, if ${\ell }_k$ is separated into real and imaginary parts as ${\ell }_k=\frac{1}{2}(a_k+\text{i}{b}_k)$, then $a_{-k}=a_k$ and $b_{-k}=-b_k$, hence Eqs. (2.1), (2.3) and (2.4) can be expressed in real-valued terms as

When the terms in Eq. (2.7) beyond $k=B$ all vanish, there are evidently $2B+1$ coefficients that can be seen to be determined exactly by using Eq. (2.9) with any $J\ge 2B+1$. (This follows more simply in the complex-valued representation where trigonometric product-to-sum rules are trivial.) In other words, the number of coefficients that can be determined exactly by way of a sum then matches the number of sample values of the original function. This sets an unbeatable efficiency target for the projections sought in the following sections. Because symmetric functions appear below, i.e. $L(-\theta )\equiv L(\theta )$, notice that only $a_k$ is non-zero in such cases and can be determined by a collapsed sum over just $0\le \theta \le \pi$. Again, the number of fitted coefficients can readily be seen to then match the number of samples.

3. Rotationally symmetric radial fit (i.e. m = 0)

To account explicitly for the role of the best-fit sphere, a pre-factor is included on the $m=0$ term in Eq. (1.1) so this component must be handled separately. It is shown in [8] that the fit required for Eq. (1.9) can be achieved in two remarkably simple steps. The first is to perform a projection in terms of an auxiliary set of orthogonal polynomials, which leads to a set of coefficients given simply by

For Eqs. (3.1) and (3.2), notice from Eq. (1.8) that

4. Radial fit for m > 0

With the estimates derived in Sec.2 for $A_m(u)$ and $B_m(u)$ of Eqs. (1.10a,b), the final step in the projection is to determine $a_n^m$ and $b_n^m$. These equations both have the same form, which can be written generically as

The determination of $c_n$ in Eq. (4.1) can be expressed in terms of a sequence of elementary operations that involve matrices with just one off-diagonal band. It is convenient therefore to define an upper-triangular matrix, say ${\text{U}}_{\text{pq}}$, by

In this notation, it is shown in the Appendix that the coefficients in Eq. (4.1) can be determined by a sequence of these simple recurrence-based operations, namely

5. Demonstration

An example of a freeform surface from the patent literature was considered in Sec. 2.2 of [6]. This surface has a plane of reflection symmetry and its conventional description involves Cartesian polynomial terms of the form $x^j{y}^k$ where $j+k\le 10$ and $j$ is always even on account of the symmetry and coordinate choice. The end result of expressing this surface in terms of $Q_n^m$ was also given in Fig. 4 of that work, but without derivation. An efficient means to achieve that conversion has been laid out above and it is now helpful for code verification to note some of the intermediate results. To give extra digits for numerical checks, all units in this section are expressed in femtometres (fm). For confirmation of the sag itself, now written with Cartesian arguments as $z=F\left(x,y\right)$, some sample values are

The symmetry in this example means that it is convenient to choose the coordinate relation $(x,y)={\rho }_{\max }(u\sin \theta ,u\cos \theta )$ so that $B_m(u)\equiv 0$, hence $b_n^m\equiv 0$. With this, for each of the azimuthal averages [such as in Eqs. (1.5) or the determination of $A_m(u)$ for Eqs. (3.4) and (4.1)] it is simplest to sample uniformly at ${\theta }_j=(j-\frac{1}{2})\pi /J$ for $j=1,2,\mathrm{...}J$. It turns out that $J=11$ is sufficient for this tenth-order surface. When $N=4$ and $K=N+2$ is used in Eqs. (4.1) and (4.8), Eq. (4.7) reproduces the coefficients in Fig. 4 of [6] to sufficient accuracy to match the sag at the sub-nm level. This involves sampling at the 66 locations shown in Fig. 1 (associated with ${\theta }_j$ and $u=\cos {\varphi }_k$ where $j=1,2,\mathrm{...}J$ and $k=1,2,\mathrm{...}K$) and yields $a_n^m$ for $m=0,1,2,\mathrm{...10}$ and $n=0,1,2,\mathrm{...4}$. Of those 55 fitted coefficients, only the 34 of order up to ten need to be retained and, as discussed in [6], these correspond to the original tenth-order Cartesian terms given in the patent. (The 21 higher-order terms that were dropped make only sub-nm contributions.)

 

Fig. 1 Sample locations for a plane-symmetric conversion with eleven spokes and six rings. Reflection symmetry means that it is sufficient to sample in just one half of the aperture.

Download Full Size | PDF

The effects of aliasing are suppressed by using larger values for $J$ and $N$ (hence also larger $K$ given $K=N+2$). Because the projection developed above is so efficient, there is little need to run with minimal values of these parameters for simple sag conversions: computing more coefficients allows confirmation of the decaying spectrum and then only the terms of significance are ultimately retained. Increasing $J$ has negligible impact in this example, but accuracy is further improved by using $N=9$ (so $K=11$). With this choice, when $m=1$, the vectors written in the Appendix as $r$, $e={\text{U}}_{\text{hk}}^{-1}{\text{L}}_{\text{hk}}^{-1}r$, $d={\text{U}}_{\text{st}}^{-1}e$, and $c={\text{U}}_{\text{fg}}d$ are found to be

The intermediate results given in Eqs. (5.2–4) are offered as an aid in the identification of coding errors. Ultimately, however, there is an acid test for verifying any code that implements such projections: plot and analyze the difference between the original sag expression and the new one. Similarly, there are analogous intermediate checks such as plotting the difference between the left and right sides of Eq. (A.9), etc. Notice too from Eqs. (5.2–4) how rapidly the coefficients decay and therefore what it means to say that such cases are “effectively band limited” (hence that aliasing can readily be suppressed). It is this property that makes it possible to use highly efficient sums in place of integrals in Eqs. (2.9), (3.5), and (4.8).

6. Fitting mid-spatial frequencies

As pointed out in Sec.3 of [5], it is sometimes helpful to express the elements within the brackets on the second line of Eq. (1.1) in terms of their amplitude and phase, say

6.1 Projections of pure sinusoids

Consider, for example, a sag function that is perfectly sinusoidal. Because the orientation of the variations influences only ${\varphi }_n^m$, it can be chosen arbitrarily and the invariant values of ${\alpha }_n^m$ are of primary interest. To be consistent with the example in Sec.5, the sag is oriented to be symmetric in $x$ and $(x,y)={\rho }_{\max }(u\sin \theta ,u\cos \theta )$. As a result, $b_n^m\equiv 0$ and it follows that ${\alpha }_n^m=\left|a_n^m\right|$ with ${\varphi }_n^m=[1-\mathrm{sgn}(a_n^m)]\pi /2$. Start then with a sag of unit peak value of the form

To clarify some of these observations, specific results are presented in Fig. 2 for the case $C=1$ where, to flatten the best-fit sphere of Eq. (1.5) and suppress the cosine effects of Eq. (1.6), I choose ${\rho }_{\max }={10}^5$. (When the sag variation is much smaller than the aperture size, the part is essentially flat, so this aperture size is used in all the simple sinusoidal cases considered here.) Notice for the anti-symmetric case, i.e. for $\delta =0$ as presented in the left column of Fig. 2, that the polynomial is dominated by the tilt, coma, and trefoil terms, i.e. $(m,n)=(1,0)$, $(1,1)$, and $(3,0)$, respectively. The symmetric case, where $\delta =\pi /2$, is dominated by spherical aberration, astigmatism of second and fourth order, and tetrafoil, i.e. $(m,n)=(0,0)$, $(2,0)$, $(2,1)$, and $(4,0)$, respectively. In both cases, the magnitude of the largest tabulated coefficients is of the order of unity and the overall fit error in both cases is below 10⁻⁶ times the peak-to-valley of the sag. The PSS is evidently almost all at third order for the odd case, and at second and fourth order for the even case.

 

Fig. 2 A sinusoidal sag function is presented in each column. Both correspond to a single full period, i.e. C = 1 in Eq. (6.2), but one is symmetric while the other is antisymmetric. The sags are plotted in the first row with a linear color ramp between the displayed extreme values. The peak-to-valley and standard deviation are also shown. The spectra of fit coefficients are tabulated in the second row as functions of m and n and include all terms up to Cartesian order t = 11. The PSS’s are presented in the third row. The rms gradient is given by the square root of the sum of these values, and the result is shown as “rss” in those plots since it is the root-sum-square of all the tabulated coefficients.

Download Full Size | PDF

When expressed in terms of a normalized coordinate, say $v$, and averaged uniformly over an integral number of periods, the rms slope of $\sin (\pi Cv+\delta )$ is $2^{-1/2}\pi C\approx 2.22C$. For integral values of $C$ therefore, although it corresponds to a differently weighted average, the number given as “rss” in the plots in the third row can be expected to be approximately equal to $2.22C$. This estimate also applies when averaging over any region much larger than a period, so the resulting approximation should be valid for all $C$ much larger than unity.

It follows from Eq. (6.2) that the results for a sinusoid of different phase can be found by simply taking a linear combination of the coefficients in the second row of Fig. 2. For example, a phase of $\delta =\pi /4$ corresponds to an equal mixture of these even and odd cases, which interleaves the results so that the entries are then all non-zero and scaled by $\cos (\pi /4)=\sin (\pi /4)=2^{-1/2}$. Similarly, the PSS is found by interleaving the two plotted partially summed spectra. Keep in mind too that rotating the sag (i.e. giving the colored contours some other orientation in Fig. 2) means that $b_n^m$ then becomes non-zero in keeping with the comments made about rotation following Eq. (6.1), but ${\alpha }_n^m$ remains untouched. For any fixed frequency $C$, therefore, the entities ${\alpha }_n^m$ and $S_t$ for the case $\delta =\pi /4$ give access to a general characterization of such a sinusoidal sag of arbitrary orientation and phase. Some results are presented in Fig. 3 for the cases $C=5$ and $C=25$.

 

Fig. 3 Sinusoidal sag functions presented much as in Fig. 2 but now as a combination of the even and odd cases and for C = 5 and 25. The fit coefficient amplitudes are plotted on a logarithmic color scale in the second row as functions of m and n. Their peak value is shown next to the red block of the legend and each colored block below that represents a reduction by a factor of ten down to 10⁻⁷ times the maximum. The third row holds the PSS’s.

Download Full Size | PDF

The spectral amplitudes, i.e. ${\alpha }_n^m$, are plotted in the second row of Fig. 3 because there are too many coefficients to tabulate. The color scale for these plots is logarithmic and ranges from red for the peak value down through seven orders of magnitude to purple. The coefficients of smaller value are not visible in this representation. In fact, ${\alpha }_n^m$ rapidly falls by orders of magnitude below the blue diagonal zone in those plots. It is striking that the contours of ${\alpha }_n^m$ coincide closely with contours of $t$ and that the highest ridge of that distribution occurs near $t\approx \pi C$. As a result, the PSS —found by summing along contours of $t$— peaks strongly near $t\approx \pi C$. More importantly, the majority of the PSS falls within a narrow interval around this peak. To emphasize this and clarify the trend, the PSS for $C=100$ with $\delta =\pi /4$ is plotted in Fig. 4. The rss values shown in Figs. 3 and 4 confirm that, as anticipated in the discussion of Fig. 2, the rms gradient is approximated well by $2.22C$.

 

Fig. 4 The PSS for C = 100 follows upon fitting over 150,000 coefficients.

Download Full Size | PDF

6.2 Frequency filtering the Q-spectra for a sinusoid

It turns out that 60% of the blue area in Fig. 4 falls within the last four lobes of the distribution, i.e. $280<t<320$, so the rss of these components alone is over 75% of the total. More generally, sinusoidal variations with a period of $\lambda \approx 2{\rho }_{\max }/C$ are principally localized in their PSS around the region $t\approx \pi C$. By combining these two relations, it follows that the abscissa for the PSS can be loosely interpreted as a spatial frequency through

A more complete appreciation of the associated frequency resolution follows upon separately summing different bands of these spectra. As an example, three sub-bands of the spectrum in the right-hand column of Fig. 3 are shown in Fig. 5. Notice that the spacing between the ridges that run from top to bottom near the centre of each plot is consistent with Eq. (6.3) given that the dominant components in each are near $t$ equal to 30, 60, and 75, respectively. It is striking, however, that spatial variations are judged by this basis to be of lower frequency when they are nearer and parallel to the edge. That is, the dominant variations near the left and right edges of these plots appear in the two lower frequency bands (which fortuitously means that a part’s typical edge effects can be more readily captured in such spectra). Most of the original sinusoid is clearly in the rightmost plot in Fig. 5, however, and the rms slope of these filtered components doubles at each of the two steps from left to right in the figure. Together, these components constitute the associated sinusoidal form to within ± 0.0017; the balance falls in $t>90$ and that residual drops by more than a factor of ten with each few additional orders in $t$.

 

Fig. 5 The data shown at top right in Fig. 3 is presented here after a low-, medium-, and high-pass filter has been applied in terms of the Cartesian order in the Q basis. These plots contain only the components in the data from (a) t = 1 to 30 at left, (b) t = 31 to 60 in the central plot, and (c) t = 61 to 90 at right. The rss (i.e. the rms gradient) for these three bands is 12.1, 24.9, and 48.5, respectively, and the standard deviation of each band appears in the plots.

Download Full Size | PDF

These Q-spectra can be filtered not just in terms of $t$, but also in $m$ (for separation of different angular frequencies) and $n$ (for different radial frequencies). An example is given in Fig. 6 where the filtering is based on angular frequency. Because the results in Fig. 6 are a combination of filtering in $t$ and $m$, they are also effectively filtering in $n$. That is, as can be appreciated by tracking the variations along either some diametral lines or concentric circles, both the radial and angular frequencies change as expected in each of the two steps from left to right. The capabilities of such processes are explored by analyzing more complex data in the next section.

 

Fig. 6 The data in the rightmost plot of Fig. 5 is now passed through a low-, medium-, and high-pass filter that discriminates based on the azimuthal order rather than the Cartesian order. These plots contain only the components in the data from (a) m = 0 to 25 at left, (b) m = 26 to 50 in the central plot, and (c) m = 51 to 90 at right. This division is chosen so that there is roughly an equal number of terms in each of these three sub-bands [two of the domains are parallelograms and one a trapezoid in the (m,n)-plane]. Together, these components exactly comprise the original function.

Download Full Size | PDF

6.3 Demonstration with synthetic data

A better appreciation of the ideas in the previous subsection can be obtained by analyzing some synthetic data that mimics as-built part shapes. For this, a superposition of a sum of thousands of sinusoidal terms with a randomized, decaying Fourier spectrum is overlaid with some additional structure. The fitting process for this case determines $a_n^m$ and $b_n^m$ over $n$ from 0 up to $N=75$ and $m$ from 0 up to $M=2N=150$. To achieve this, samples are taken at $u=\cos {\varphi }_k$ with $k=1,2,\mathrm{...}K$ where ${\varphi }_k$ is given in Eq. (3.6) and $K=N+2=77$. Since there is no plane of symmetry, the sampling is no longer in a half-plane as in Fig. 1. Instead, the azimuthal averages used to determine $A_m(u)$ and $B_m(u)$ for Eq. (4.1) involve sampling at $\theta =j\Delta$ for $j=1,2,\mathrm{...}J$ where $\Delta =2\pi /J$ and $J=2(M+1)=302$. That is, 302 samples are taken on each of 77 rings so, in this case, $77\times 302=23,254$ sag samples are used to determine $(N+1)\times (2M+1)=22,876$ coefficients. This is close to the one-for-one efficiency standard set by the discrete Fourier methods. As an aside, since the Fourier series on each ring is independent, an attractive and generally equivalent alternative is to sample on the $k\text{'}\text{th}$ ring at $\theta =(j+k/\phi )\Delta$ where $\phi =(5^{1/2}+1)/2\approx 1.618$ is the golden ratio. The distinction is displayed in Fig. 7 and is related to Fibonacci lattices [11]. Of course, the number of samples around these rings can be boosted to a power of two to increase the benefit delivered by FFT algorithms, if desired [10].

 

Fig. 7 The default sample locations at left for N = 25 (hence M = 50, J = 102, and K = 27) involves 2,754 samples to extract (2M + 1)(N + 1) = 2,626 coefficients. This is appropriate for a general fit of a spectrum like that in the left column of Fig. 3. The Fibonacci-inspired and analytically equivalent option at right avoids unattractive bunching in the outer zones.

Download Full Size | PDF

The spectrum for this example is plotted in Fig. 8 and the data itself for the normal departure appears in Fig. 9(a). At a demanding yet realistic level, if the units are taken to be millimeters and the departure needs to be characterized down to the nanometer level, the fit must be accurate to better than one part in 10⁶. In this case, it is evident in the plot at left in Fig. 8 that only about a half of the 23,000 fitted coefficients are large enough to be significant for these purposes, viz. those within $t\le 140$. (Remember that these fits must be taken to sufficiently high orders to suppress corruption by aliasing, thus subsequently disregarding higher orders is an essential aspect.)

 

Fig. 8 The peak coefficient in this spectrum is 2.3 and the color legend steps down to purple for an amplitude of 2.3e-7. The PSS is also plotted on a log scale as the green curve at right, while the blue curves give the results when the sampled sag values are corrupted by additive noise uniformly distributed between ± 1.0e-3, ± 1.0e-4, and ± 1.0e-5. Only the highest level of noise can be distinguished from the noise-free case in the left half of the plot (t < 80). The purple curve is found by using interpolation from a 250x250 pixelated grid.

Download Full Size | PDF

 

Fig. 9 (a) Normal departure extracted from a synthetic data set. The components in that data from (b) t = 1 to 10, (c) t = 11 to 24, (d) t = 25 to 64, (e) t = 65 to 84, and (f) t = 85 to 140.

Download Full Size | PDF

It is instructive now to perform some filtering. First, it is clear from the plot at left in Fig. 8 that the spectrum is dominated by the terms satisfying $t\le 10$. Figure 9(b) reveals how well the shape is captured by those first $4+2\times 30=64$ fit coefficients. The residual is almost indistinguishable from what is plotted in Fig. 9(c), which is the sum of the components in the spectrum within $11\le t\le 24$. Striking structures become visible at the higher spatial frequencies presented in Figs. 9(d), 9(e), and 9(f). Such MSF patterns are oftentimes found on parts manufactured with modern sub-aperture tools: spoking is evident in Fig. 9(d) and two sets of raster marks at 60° to each other are clear in Fig. 9(e) while Fig. 9(f) suggests some high-order rings on the surface. As a further demonstration, the spoke and high-frequency ring components are more specifically isolated in Fig. 10, but the more impressive and central result is the separation of spatial frequencies demonstrated in Fig. 9: the capability to decompose the data in Fig. 9(a) into the five bands shown is both striking and useful.

 

Fig. 10 Plots of the components of the spectrum in Fig. 8 for the data in Fig. 9(a) with (a) m = 14, 28, 42,… 140 and (b) m = 0 and n > 7. Notice that the multiples of 14 stand out as vertical bands in the plot at left in Fig. 8 and that there are anomalously large values for large n when m = 0. These are the indicators that suggest these two final filtering steps.

Download Full Size | PDF

More realistically, metrology is corrupted by noise of various types. An indication of the impact that noise has on these fitting processes is presented in the PSS in Fig. 8. It turns out that the spectral coefficients are all impacted at about the same absolute level. Further, it is pleasing that summing the difference between the clean and noisy spectra for coefficients above the noise floor yields a map whose peak value is comparable to the noise level. Another type of noise is encountered when interpolating from pixelated metrology to obtain the desired sample values. To explore this, I used a $250\times 250$ pixel grid with bi-linear interpolation and the resulting PSS is shown as the purple curve at right in Fig. 8. In this case, summing all terms within $t\le 140$ gave an absolute error less than 0.0005. This means that all the sub-bands in Fig. 9 except (f) can be extracted in this way. Interestingly, Fibonacci sampling yielded equivalent results, but gave slightly larger coefficients at high azimuthal orders for this interpolated case.

For parts with circular apertures, sampling processes of the type shown in Fig. 7 appear to be promising partners for coordinate measuring machines. When the metrology yields a pixelated map on the other hand, although interpolation is then a workable option, the process can in practice be complicated by the presence of regions of invalid data. In such cases and, more generally, for intentionally non-circular apertures there are useful modifications based on least-squares methods for the fitting processes developed above. While the factoring into azimuthal and radial stages remains a critical step, the details are not discussed further here.

7. Concluding remarks

The goal of this work was to develop a robust projection for expressing nominal shapes in terms of the gradient orthogonal basis $Q_n^m$. Although the natural projection would involve derivatives for determining gradients, this could sometimes be problematic and was found to be unnecessary. Despite the fact that the underlying orthogonality is based on an integral, the primary results, namely Eqs. (3.3), (3.5), (4.7) and (4.8), are of comparable efficiency to discrete Fourier transforms. That is, they yield almost one fitted coefficient for every sample point evaluated within the circular aperture. Further, the first stage of this projection can be implemented by using FFTs, and the second uses Chebyshev-Gauss quadrature to deliver its remarkable overall efficiency. It is also worth noting that the projection developed in the Appendix to an intermediate basis can equally be followed by a change of basis using Eqs. (B.7-9) of [5] together with Salzer’s recurrence-based methods of [12]. In this way, the results can be expressed in terms of other orthogonal polynomials if required, e.g. Zernike polynomials (although there can be stability problems at extreme orders in that case). That is, the methods developed above allow freeform shapes that are defined over a disc to be expressed in terms of a variety of standard orthogonal polynomials.

Because they readily extend to many thousands of terms, the methods developed here create possibilities for a novel characterization of MSF on as-built parts —regardless of whether the parts are nominally rotationally symmetric or intentionally freeform. This has long been thought to be beyond the reach of polynomials alone. It is also significant that early innovations for coping with the evolving challenges of tolerancing MSF led to a realization of the importance of an rms gradient metric in some applications [13]. It is noteworthy then that the projection to $Q_n^m$ gives direct access to the rms gradient, and that this is connected to the rms radius of a ray spot diagram. Further, the power of the spatial frequency-filtering processes demonstrated in Sec. 6 is perhaps the most promising surprise to emerge from this work. It establishes that, on top of their benefits in nominal shape characterization and design, tailored orthogonal polynomials also have the potential to serve as a useful diagnostic and qualification tool for MSF-related aspects of optics production.