Characterizing the shape of freeform optics

: A recently introduced method for characterizing the shape of rotationally symmetric aspheres is generalized here for application to a wide class of freeform optics. New sets of orthogonal polynomials are introduced along with robust and efficient algorithms for computing the surface shape as well as its derivatives of any order. By construction, the associated characterization offers a rough interpretation of shape at a glance and it facilitates a range of estimates of manufacturability.


Introduction
Progressive ophthalmic lenses were conceived over 75 years ago [1] and soon became more complex [2,3]. They clearly demonstrated the unique benefits that freeform optical surfaces can offer through the extra degrees of freedom available in surface shape. As documented in [4], other prominent commercial applications emerged in the 1970's. On account of advances in optical fabrication and metrology, interest in a wider class of applications has escalated in recent decades. Within the domain of precision optics, these applications include headmounted and head-up displays, phase elements for computational imaging and unobstructed reflective imaging systems (from ultra-compact digital cameras through projection television to EUV lithography). Lower-precision freeform surfaces have long been prominent in nonimaging applications such as illumination and it is impressive that some of the associated design methods have expanded in their capabilities so that they are now valuable for imaging applications as well [5,6]. Other interesting design methods are also evolving to enable the discovery of more effective systems involving freeform optics [7][8][9].
Standardized methods for communication of the nominal shape of freeforms are critical to this technology. Although no single convention can meet all the requirements of the various contexts and applications, some reasonably general-purpose options exist, e.g. point clouds with interpolation and various splines or wavelets [e.g. 10,11]. Alternatively, a common method in precision optics is to express the surface sag explicitly as a conic base with an additive deviation expressed as a polynomial. This has the advantage of being infinitely differentiable and is a simple generalization of the most widely used convention for characterizing rotationally symmetric optics. When the polynomial basis is orthogonalized, the associated coefficients become a spectrum of the surface's frequency content that can be interpreted readily upon inspection. Further, such a spectrum can also facilitate estimates of manufacturability (which typically involve consideration of generating, testing, and oftentimes polishing). When implemented by using recurrence relations, this option is so robust and efficient that it enables an unlimited number of coefficients to be used as needed. These observations motivated the development of the polynomials introduced in this work.
Basic geometric ideas are presented in Section 2 where a general form is proposed to characterize the sag of a freeform optic. The polynomials used in this sag expression are derived in Section 3 and algorithms for their robust computation are presented in Appendices A and B. For demonstration, a particular surface is expressed in these terms in Section 4. A generalization is offered for certain special cases in Section 5 where some simple measures of testability are also discussed before the concluding remarks of Section 6.

Coordinates, aperture shapes, and sag expressions
Off-axis sections of rotationally symmetric surfaces can be regarded as a first step towards freeform optics. In such cases, it is natural to adopt coordinates that are aligned with the part's axis of revolution. For example, the sag of a conventional conic section of revolution can then be expressed as z = c con [(x + s) 2 + y 2 ] 1+ 1− (1+κ )c con 2 [(x + s) 2 + y 2 ] = c con (ρ 2 + 2sρ cosθ + s 2 ) 1+ 1− (1+κ )c con 2 (ρ 2 + 2sρ cosθ + s 2 ) , (2.1) where c con is the conic's axial curvature and κ the conic constant. The first expression in Eq. (2.1) involves the Cartesian coordinates {x, y, z} with their origin displaced from the part's vertex by s ; the second is in terms of the associated cylindrical polar coordinates {ρ,θ, z} where x = ρ cosθ and y = ρ sinθ . Such a configuration with a vertical z axis is sketched in Fig. 1.
For freeform optics that are to be turned/ground, smoothed, polished, and tested -perhaps interferometrically with nulls or near-nulls and/or stitching-the more cost-effective shapes tend to be closer to spherical. Manufacturability typically grows in difficulty with departure from a best-fit sphere. More specifically, it is the rate of change of the departure along the normal from the best-fit sphere that must often be considered in both fabrication and testing. As discussed in [12], the difference between the local principal curvatures on the surface is often a driving consideration that can be related directly to derivatives of the normal departure. For many purposes, it is therefore helpful to characterize freeform shape by explicit reference to a best-fit sphere. That is one of the key motivators for the developments reported in what follows. Another of the critical considerations is that, for efficiency and numerical robustness, an orthogonal basis delivers significant benefits when characterizing shape [13,14]. Fig. 1. An off-axis section of a paraboloid is commonly encountered and has some of the challenges of more general freeform optics. The coordinate origin used for Eq. (2.1) is neither at the top nor the bottom of the cylinder as drawn for clarity at left, but in a plane that is tangent to the part's vertex as shown at right. The separation between the cylinder and part axes is s.
It is not only the surface's shape but also the shape of the part's perimeter that must be characterized to fully specify a freeform. Rectangles, hexagons, ellipses, and kidney beans are familiar examples of nominal aperture shapes, but their specification is not always as simple as it may seem because the part's perimeter typically does not lie in a plane. The aperture shown as the red curve in Fig. 1, for example, is defined as the intersection of a circular cylinder with the parent surface. For other geometries, the proposal in this work is to tightly enclose the part within a circular cylinder and then use the polar coordinates associated with that cylinder in the characterization of the surface's shape, see Fig. 2. Orthogonalization is necessarily performed over a specific domain and, for that purpose, I adopt this circular region that encloses the part. Fig. 2. The characterization of a segment of a general freeform surface is depicted at left. In this case, as may be familiar from segmented mirrors, the segment of interest is taken to be the intersection of a hexagonal column with the surface. I opt for the coordinate origin to sit on the surface and choose the green cylinder to tightly enclose the segment of interest on the surface. The cylinder will typically be chosen to be nominally normal to the surface segment, as shown. The case for a part that resembles the paraboloidal segment of Fig. 1 is depicted at right. The next step therefore is to extend Eq. (2.1) of [14]. In particular, if we define u by u := ρ / ρ max where ρ max is the radius of the tightly enclosing cylinder, then the sag is now Here, Q n m (v) is a polynomial of order n in v , c is the curvature of what I somewhat loosely refer to as the "best-fit sphere", and 1− c 2 ρ 2 is the cosine of the angle between the z axis and the normal to that sphere. The presence of the cosine factor in Eq. (2.2) means that, to a first-order approximation, the entity within braces corresponds to a departure from the sphere that is measured along the local normals to that sphere. Notice that the sum on the first line of Eq. (2.2) has precisely the same form as that used in [14] except that the coefficients and polynomials now carry a superscripted zero. This zero distinguishes them from the non-rotationally symmetric components, namely a n m , b n m , and Q n m (u 2 ) with m > 0 . Also note that the average over θ of f (ρ max ,θ ) is precisely the sag of the best-fit sphere at ρ = ρ max : the component inside the braces averages to zero. In the unusual event of having to fit a predetermined shape with this representation, this observation can be useful to evaluate c. More typically, all of the coefficients in Eq. (2.2) will be determined through system optimization during design.
Although the form of the polynomials is not locked down until Section 3, the similarity between the sum on the second line of Eq. (2.2) and the familiar decomposition into Zernike polynomials establishes that the terms in braces provide a complete set. As can be shown by using trigonometric identities, these terms are just a particular combination of the Cartesian elements X j Y k where X = ucosθ and Y = usinθ . [It is useful in this to equate real and imaginary parts after expanding the identity u m (cos mθ + isin mθ ) = ( X + iY ) m .] The maximum order of the elements that compose the contribution associated with a n m or b n m where m > 0 is just m + 2n, i.e. this is the maximum power of u in those contributions. It would also be natural, therefore, to truncate the sums on the second line of Eq. (2.2) to include only the terms for which m + 2n ≤ T for some T . Due to the prefactor on the sum in the first line, where m = 0 , it is consistent for that sum to then extend up to 2n + 4 ≤ T .
Without loss of generality, it is possible to assume in Eq. (2.2) that When they're non-zero, the contribution from these terms to the normal departure is proportional to u (a 0 1 cosθ + b 0 1 sinθ ) ∝ (a 0 1 x + b 0 1 y) . It is clear therefore that the cylinder's axis can be chosen to remove these average-tilt terms (much as the best-fit sphere was chosen to make the rotationally symmetric component vanish at the cylinder's centre and edge). Alternatively, the cylinder can be oriented to remove either the mean tilt around its rim or the local tilt at the origin, which lead respectively to a n 1 Q n Neither of Eqs. (2.4a) and (2.4b) is functionally as convenient as Eq. (2.3). I mention these options mainly to highlight some of the residual freedom in configuring the cylinder. Any of them can be used during optimization as part of design provided that the surface has appropriate tip and/or tilt parameters as separate degrees of freedom. In this way, it is assured that the base sphere of Eq. (2.2) is a decent approximation to the best-fit sphere as far as manufacturability is concerned. In the terminology of metrology, for example, the irrelevant piston, tilts, and power are suppressed in this way. The most important issue that remains is the process of constructing the new polynomials, i.e. Q n m (v) .

Orthogonalization in terms of the mean square gradient
Because it is the rate of change of the normal departure from a best-fit sphere that often drives manufacturability, the polynomials discussed in [12] were constructed to be orthogonal in slope. The mean square slope is then given by the sum of the squares of the coefficients in the departure. For freeform optics, it is the modulus of the gradient of the normal departure from the best-fit sphere that plays this role. That is, if the term in braces in Eq. (2.2) is written as the natural extension of the basis used in [12] is to construct Q n m (v) so that the mean square gradient of this normal departure from the best-fit sphere satisfies ( ) Here, b n 0 = 0 and angle brackets denote the average over the aperture defined by The weight function used in [12] Equations ( 2) leads to sums of products, and the product of any two distinct azimuthal orders averages to zero. The only remaining step is to choose the polynomials within each separate azimuthal order so that they are internally orthonormal. Of course, Q n 0 (v) are precisely the ("Qbfs") polynomials introduced in [13], but those for m > 0 require further treatment. After performing the angle integral, the condition for orthonormality reduces to the requirement that is unity when n = ′ n and vanishes otherwise. With this, it is straightforward to construct the polynomials by using Gram-Schmidt orthogonalization. The results for m from 0 to 5 and n from 0 to 3 or 4 are presented in the tableau in Fig. 3. Plots of the cosine versions of the associated functions in Eq. (3.1) along with contour plots of the modulus of their gradient are given in Fig. 4. The complementary basis members, i.e. u m sin(mθ )Q n m (u 2 ) , have similar plots, of course, but each is rotated by π / (2m) .  Fig. 1 of [14], although now it is the absolute value of the radial slope that is plotted. Recall that it is the non-uniform weight used in Eq. (3.4) that helps to constrain the peaks in the magnitude of the gradient to be not much bigger than 2 . This peak value is in keeping with a sine-like shape that has an rms value of unity. Also keep in mind that, as shown in Fig. 5, the gradient is a vector and it is the integral of the inner product of these vector fields that underpins the orthonormalization. The polynomials listed in Fig. 3 allow for over 50 coefficients to be used in characterizing a freeform's shape [although the truncation discussed in the paragraph before Eq. (2.3) means that a subset of only 1+ 2 × (3+ 2 + 2 + 1+ 1) = 19 of these would be used if we stop at T = 5, and Eqs. To go beyond low orders, recurrence relations offer efficient, robust algorithms for evaluating arbitrary linear combinations of Q n m and of their derivatives. One option to achieve this involves a link to Jacobi polynomials, and the associated details are presented in the Appendices. It is worth noting that a similar set of polynomials is developed in [15] to provide a basis for curl-free vector functions by way of a set of gradient-orthogonal polynomials.
Interestingly, it is the different weight function [i.e. w(u) in Eq. (3. 3)] adopted here that causes the main difference with the polynomials derived in that work. It is sometimes helpful to notice that the component in square brackets in Eq. (3.1) can trivially be re-expressed in terms of an alternative set of coefficients that give amplitude and phase: a n m cos(mθ ) + b n m sin(mθ ) = α n m cos(mθ − φ n m ). (3.6) With this, the mean square gradient of Eq. (3.2) satisfies In this sense, α n m is then of primary importance in the spectrum while the phases become secondary. That is, a rough assessment of the dominant components of a part's shape then follows upon inspection of just one half of the original number of coefficients.  Fig. 4. It is the point-by-point inner products of different members of such maps that integrate to zero. The scale is fixed in all of these plots so that an arrow whose length matches the grid spacing represents a gradient of magnitude equal to the square root of two.

Simple surface for demonstration
For demonstration and code verification, a simple paraboloid is expressed in the terms proposed above. This was also done in Section 4 of [14], but now an off-axis section of that same paraboloid is fitted with respect to a locally best-fit sphere. The paraboloid has curvature 1/ c parab = 20mm and the distance from its axis to its point of intersection with the centre line of the green cylinder drawn in Fig. 7 is taken to be s = 20mm . When the radius of the cylinder is ρ max = 10mm and the cylinder axis is chosen to be normal to the surface at its point of intersection, the curvature of the best-fit sphere is given by 1/ c = 37.405283mm . (This curvature is given to an artificially high precision for code verification of a nm-level fit of the paraboloid with respect to this sphere.) The fit results are shown in Figs. 6 and 7. Notice that symmetry in y means that b n m ≡ 0 . These results use some additional terms beyond those given in Fig. 3: While the coefficients in Fig. 6 correspond to a configuration satisfying Eq. (2.4b), the results for the configuration satisfying Eq. (2.3) are given in Fig. 8. Notice the similarity between the two sets and that the shape is largely dominated by about 600µm of astigmatism ( m = 2 and n = 0 ) and 200µm of coma ( m = 1 and n = 1), see Fig. 4. That is, the character of the departure plot in Fig. 7 can be anticipated by a glance at these tables of coefficients. Fig. 7. The sample paraboloid is shown at bottom along with the cylinder that encloses the offaxis section of interest and the associated best-fit sphere in red. The sag departure is plotted in the upper row along with the difference between this sag and that associated with the coefficients of Fig. 6. Notice that the fit error remains less than a nanometer. Fig. 8. The corresponding coefficients for the case of Fig. 6, but where the cylinder has been tilted by 20.2223 mrad to enforce Eq. (2.3). That is, as highlighted above in red, a 0 1 now vanishes. The inverse of the best-fit curvature for this case is 37.432729 mm.

5.2)
To within a first-order approximation, the additive departure in sag from the conic in Eq. (5.1) can be converted to the departure along the normal to that conic by multiplying it by the cosine factor σ (ρ,θ ) . It follows that δ (u,θ ) itself can again be regarded as the departure along the normal to the reference surface. This generalization enables the handling of extremely fast cases that would otherwise involve unworkable hyper-hemispheres. All the methods derived in the appendices carry over without change because, once again, it is the sum in Eq. (B.2) that stands as the core component. This is also an effective option when characterizing parts that are close in shape to a conic section for which an interferometric null test is in hand and where the metrology represents a critical part of manufacturability. A part that approximates a paraboloid, for example, can be tested by using a retro sphere and a transmission flat as depicted in Fig. 1 of [12], but now imagine that only the upper half of that figure holds the off-axis segment of interest. As described in [12], a map of the fringe density in the resulting interferogram is approximately proportional to ∇δ (u,θ ) times the cosine of the angle of incidence of the test rays at the part. In many cases, the cosine factor is roughly constant. While it is the peak fringe density that ultimately drives testability, the rms fringe density is a useful and more readily accessible measure. Thankfully, the weight function of Eq. (3.4) serves to constrain the ratio of peak to rms values. When such a test is performed at wavelength λ on an N × N pixel grid, it follows from the discussion in [12] that the rms fringe density in Nyquist units, say γ , satisfies Results such as this can be used to derive simple design constraints and valuable estimates of testability.
It is striking that Eq. (5.3) involves little more than summing the squares of the shape coefficients; there is no need to evaluate a single departure from a best-fit conic or to average its squared rate of change. The orthogonalization introduced in Section 3 was constructed to make this rms gradient trivially accessible. A similar result applies, of course, for a conventional non-null interferometric test with a transmission sphere for a part described by using Eq. (2.2). [The single bounce off the part in that case means that the factor of 16 in Eq. (5.3) is replaced by an 8.] What is more, the same analysis also applies when considering the line spacing of a CGH null. Do not forget, however, that the sum of the squares in Eq. (5.3) gives an rms over the enclosing cylinder rather than just the part's clear aperture. In the later stages of design or fabrication, more accurate analyses will typically be required.

Concluding remarks
The diversity of applications of freeform optics means that no single convention for characterizing shape will be ideal for all cases. For example, when a part has a rectangular aperture and is manufactured by rastering processes, it would be natural to construct a description in terms of Cartesian coordinates rather than the cylindrical polar coordinates used above. Alternatively, for parts that are ready to be deployed immediately after generation by single-point diamond turning, it is the azimuthal component of the gradient of the sag that oftentimes drives manufacturability, and direct access to this property could be factored into their characterization. However, for high-precision optics that must be turned/ground, smoothed, finished, and perhaps interferometrically tested (probably with nulls or near-nulls and/or stitching), it is effective to express their shape with respect to a best-fit sphere. In particular, it is typically the gradient of the normal departure from the best-fit sphere that is of crucial importance. This observation motivated the generalization of Qbfs (of [13]) that was presented in Sections 2 and 3 in terms of an enclosing cylinder. The further generalizations offered in Section 5 include not only a conic constant, but also the option to move the base surface's center of curvature away from the cylinder's axis. This second freedom can be exploited even when κ = 0 . For general purposes, however, it is Eq. (2.2) that is preferred over Eq. (5.1) because other aspects of manufacturability are more readily assessed when a freeform surface is expressed with respect to such a best-fit sphere.
When supplemented with the powerful algorithms presented in the appendices, Eq. (2.2) offers a general-purpose option for characterizing freeform shape in terms of coefficients that form an intuitive and elegant spectrum. For ray tracing through such a freeform, it can be helpful to have explicit access to the surface's local normal vector. In terms of coordinates related by {x, y, z} ≡ {ρ cosθ, ρ sinθ, z}, the Cartesian components of the normal vector for the surface described by z = f (ρ,θ ) are readily found to be where subscripts denote partial derivatives. Similarly, for first-order aberration analysis, it can be helpful to note that this surface can be approximated near the coordinate origin by The results presented in the appendices give simple access to these linear combinations of Q n m (0) for fixed m > 0 , see Eqs. (A.23) and (B.2) to (B.4). Notice that this approach is a direct analog of Eq. (3.16) of [14], which can now be re-used for the sum with m = 0 in Eq. (6.2). This further highlights the important observation that Eq. (2.2) includes as a special case the option for characterizing rotationally symmetric aspheres that was introduced in [13].
Because it is built on notions of manufacturability, this single framework can therefore meet many of both the traditional as well as the emerging needs of our industry.

Appendix A: Auxiliary polynomials
Be warned at the outset that the equations in this opening paragraph will be modified below in order to avoid an isolated and surprising failure. The key here is that the polynomials defined in Section 3 can be expressed in terms of the Jacobi polynomials. In particular, consider the auxiliary polynomials defined by where the constants are given explicitly by Just as the auxiliary basis did for the polynomials used in [14], P n m (x) and Eq. (A.2) underpin a number of efficient processes for working with Q n m (x) . The scaling factor in Eq. (A.1) is chosen so that the peak absolute value of u m P n m (u 2 ) over 0 ≤ u ≤ 1 never exceeds about unity for all m and n. As in [14], this makes it straightforward to determine when terms are negligible and can be dropped from a linear combination.
In this case, the first two polynomials are This patch means that, when m = 1, Eq. (A.2) is valid only for n > 2 , and the initialization for that case is then If the integral in Eq.
This process involves replacing whichever of the expressions on the left-hand sides of Eqs. (A.8)-(A.12) appears first by the associated right-hand side, and the final step is to invoke Eq. (A.10) once more but with m replaced throughout by m + 1. While this sounds somewhat daunting, present-day computer algebra is up to the task when provided some guidance. In this way, it follows from the orthogonality of the Jacobi polynomials that P n m , P ′ n m is zero whenever n − ′ n > 1. The diagonal elements, say F n m := P n m , P n m , turn out to be given by n > 0 and m = 1, n > 0 and m > 1, where χ = m + n − 2 as a shorthand, δ m n is the Kronecker delta, and Similarly, the off-diagonal elements, say G n m := P n m , P n+1 m , are found to be n > 0 and m = 1, derivatives of any order were presented in [14]. For any fixed m ≥ 1, this appendix addresses the additional core task of working with the sums of the form of those on the second line, say where c n may be either a n m or b n m . (To reduce clutter, neither S nor c n takes a superscripted m.) By following the same procedure discussed in Section 3 of [14], it follows from Eq. (A.17) above that S(x) can also be expressed in terms of the auxiliary polynomials as One of the most striking aspects of Eqs. (B.6) and (B.9) is that, as described in Appendix B of [16], they allow derivatives of any order to be evaluated with similar ease. If S ( j ) (x) denotes the j 'th derivative of S(x) then it follows from Eqs. (B.9) and (B.6) respectively that S ( j ) (x) =