## 1.

## Introduction

With the advancement of modern precision optical fabrication and measurement technology, optical components with different surface types can be gradually realized. Under this condition, a spherical surface is not the only choice in current optical design. Owing to their abundant degrees of freedom and strong ability of aberration correction, optical aspheric and freeform surfaces are strong candidates. From the geometrical viewpoint,^{1} an optical freeform surface has nonrotationally symmetric features. From the aspect of fabrication and design,^{2} an optical freeform surface is regarded as an optical surface that leverages a third independent axis during the fabrication process to form the optical surface with as-designed nonsymmetric features. Both viewpoints reveal the basic feature of the nonrotational symmetry of optical freeform surfaces, which is a challenge for the optical community. In the last 10 years, the technologies of multiaxis slow-servo single-point diamond turning and computer-controlled small-lap polishing have been further enhanced, and the manufacturing capacity and metrology of nonrotationally symmetric optical surface have been improved. Therefore, the application fields of nonrotationally symmetric optical surfaces have increased considerably.

Generally, freeform optics can be categorized into nonimaging and imaging freeform optics. Applications of optical freeform surfaces in the nonimaging category include beam shaping,^{3} concentrators,^{4} and illumination,^{5} where optical freeform surfaces can control the light direction to improve energy uniformity and efficiency. Examples of applications in the imaging category are ultrashort projection lenses,^{6} head-mounted displays,^{7} ultraviolet lithographic objectives,^{8} and off-axis reflective infrared imaging systems,^{9} where optical freeform surfaces are used to improve imaging performance and compactness. Hence, freeform optics is regarded as the next-generation and ground-breaking technology in modern optics.^{2}^{,}^{10}

However, the question of how to represent or characterize an optical freeform surface remains. Because the optical freeform surface representation technique has a close relationship with its design, manufacturing, testing, and final application, it encompasses all the processes used in freeform optics. Namely, the freeform surface representation technique is a fundamental and key research topic, and its enhancement can further facilitate the development of freeform optics. Obviously, the optical freeform surface representation technique employed in different applications is not unique, and its selection depends on practical considerations. So far, XY-type, Zernike, Q-type, Chebyshev, and Legendre polynomials, the radial basis function, the nonuniform rational basis spline (NURBS) function, and other hybrid methods and techniques have been used to represent the optical freeform surface in each specific application. Thompson and Rolland^{2} described the history and the revolutionary character of optical freeform surfaces. Fang et al.^{11} presented a comprehensive review on the manufacturing and measurement of freeform optics, which contained some information on freeform surface representation techniques. Gross et al.^{12} published an overview of surface representations for freeform surfaces particularly from the aspect of power spectral density in the midspatial frequency for manufacturing. Fähnle et al.^{13} presented a special section on freeform optics to discuss the current developments in all aspects of freeform optics. We believe that an elaborate analysis of optical freeform surface representation techniques is necessary, as it can provide a large map about the development of these methods.

In this paper, we present a comprehensive review of the different types of freeform surface representation techniques and their applications. Section 2 describes several typical applications using different types of optical freeform surfaces. Section 3 discusses various mathematical functions and techniques that can be applied to characterize optical freeform surfaces and analyzes methods used to obtain a freeform surface from discrete data points. In Sec. 4, the differences between freeform surface representation techniques are discussed and compared in detail. Section 5 concludes with the state of the art in optical freeform surface representation and its future trends.

## 2.

## Applications of Different Types of Optical Freeform Surface Characterization

An optical freeform surface has the advantages of improving the imaging performance, reducing the number and weight, and further increasing the compactness of optical components, especially in off-axis imaging optical systems. This is owing to its multiple degrees of freedom and off-axis aberration correction ability. One of the most typical examples of a freeform surface was used by Maitenaz^{14} and Kanolt^{15} in progressive multifocal ophthalmic lenses, which date back to 1954. It was applied for the correction of presbyopia and offered the wearers comfortable, relatively clear vision with no separation between long and short distances,^{16} as shown in Fig. 1. The profile of the progressive lens is obviously of the nonrotationally symmetric type, which can be represented by the combination of Zernike polynomials or the B-spline function.

Another example of freeform surfaces used in an off-axis catadioptric Polaroid SX-70 camera was introduced by Plummer;^{17} this might be the first commercial product employing optical freeform surfaces characterized by XY polynomials to improve its imaging performance, as displayed in Fig. 2. Owing to the optical path configuration of off-axis refractive and reflective types, it had more compactness. However, off-axis aberrations aggravated its imaging quality. Thus, two freeform optical components were used to correct the off-axis aberration. A freeform corrector (D) with a fourth-order XY polynomial profile was located at the system stop, which could reduce the astigmatism, coma, and spherical aberration. Another freeform eyepiece lens (G) with a sixth-order XY polynomial profile was positioned at the viewfinder to control the astigmatism and field curvature.

In the last decade, with the development of advanced fabrication equipment, as well as commercial optical design tools and optical testing technology, the advantages of optical freeform surfaces have been extensively utilized in different applications such as projection optical systems, head-worn displays, and off-axis reflective optical systems. In these applications, the used freeform surface is generally a continuous and smooth surface characterized by different analytical functions. Additionally, a microfreeform-prism-based lens array is employed in the artificial compound eye^{18}19.^{–}^{20} and can realize a large field of view using multiple imaging channels with lower distortion. Each imaging channel is a microfreeform-prism lens represented by XY polynomials or Zernike polynomials to correct the aberration in the corresponding channel. This type of freeform surface is stepped and noncontinuous globally. From the perspective of the overall continuity of surface, the optical freeform surface can be categorized into two types, as described above. Continuous and smooth freeform surfaces are much more common in the design and application of freeform optical systems. Therefore, in the following, we discuss applications that use different types of optical freeform surfaces by focusing on continuous and smooth freeform surfaces.

A traditional projection system always has a relatively long throw distance, which involves large unused spaces. On the other hand, if the presenter stands between the screen and the projector, his shadow is shown on the screen. As the accessory of the traditional projection system, freeform optical elements can obtain an ultrashort throw distance, a large screen, and high performance projection, as shown in Fig. 3. Hitachi^{21} reported that a projector (HCP-A8) using a freeform reflective mirror could realize a 60-in. projection screen at a throw distance of only 47 cm; however, the detailed function type of its profile was not revealed. The simultaneous multiple surface (SMS) method^{22} was applied by Miñano and Benítez in the design of a freeform accessory optical system for ultrashort throw distance projection. The discrete data point cloud characterizing the freeform surface was obtained by the SMS method and could be reconstructed accurately by analytical polynomials or functions. Yang et al.^{23} used XY polynomials to obtain a miniaturized and short projector. Zhuang et al.^{6} applied odd XY polynomials representing the freeform reflector together with a field curvature correction method to design an ultrashort projector. The introduction of optical freeform surfaces to projection optical systems could further enhance projection display technology to achieve an ultrashort throw distance, larger display screen, higher imaging projection performance, miniaturization, and low weight.

Head-worn displays, which can interact with the virtual digital world and connect with the real world, are widely used in modern education, medical treatment, entertainment, and military training. The use of a freeform prism as the eyepiece lens for head-worn displays, as shown in Fig. 4, was reported by Okuyama^{24} and Takahashi.^{25} Each profile of the freeform prism can be represented by XY polynomials or Zernike polynomials. Cakmakci and Rolland^{26} designed a dual-element off-axis near-eye optical magnifier by combining a freeform reflector represented by XY polynomials and a diffractive lens. Furthermore, Cakmakci et al.^{27} applied the radial basis function, instead of XY polynomials, to characterize the freeform reflector, which could increase the exit pupil diameter from 8 to 12 mm. Cheng et al.^{7} reported a freeform prism-based head-worn display characterized by XY polynomials, which had the advantages of a fast focal ratio, large field of view, and see-through properties. Zheng et al.^{28} designed and analyzed an off-axis see-through head-worn display with XY polynomials that had a large eye relief of over 60 mm. Pan et al.^{29} demonstrated an off-axis two-freeform reflector described by XY polynomials with an 8-mm exit pupil diameter, which was used to enhance the visual sensing of slightly visually impaired patient. The development of optical freeform surfaces advances head-worn displays in the direction of a large exit pupil diameter, long eye relief, large field of view, high-resolution display, more comfort, and miniaturization.

Optical freeform surfaces are frequently used in off-axis reflective optical systems. Xu et al.^{30} presented a freeform spectrometer using a toroid reflective surface to correct astigmatism. Rodgers^{31} designed an off-axis four-mirror system that used the off-axis section of an even-order asphere. Furthermore, Nakano and Tamagawa^{32} designed an off-axis three-mirror system. The original configuration was analyzed by conic surfaces and further optimized by Zernike polynomials. Fuerschbach et al.^{9} presented an off-axis three-mirror long-wavelength infrared imaging system using the off-axis part of a fringe Zernike polynomial surface by combining the nodal aberration theory and the full-field display method, as shown in Fig. 5. The use of an off-axis optical path configuration and a freeform surface makes the imaging system more compact and more efficient in correcting the off-axis aberration.

Optical freeform surfaces also play a key role in panoramic imaging systems with a wide field of view. Ma et al.^{33} presented a panoramic annular imaging system with nonsymmetric, rotational, and variable focal length, which had a freeform surface characterized by XY polynomials, as displayed in Fig. 6. Zhou and Bai^{34} designed a small-distortion panoramic annular lens with a Q-type freeform surface. Previous panoramic imaging systems had an inherent shortcoming, namely, a central blindness without the image. By applying a stitched freeform surface between the central and peripheral regions, a large field-of-view panoramic lens without central blindness could be realized.^{35}

Additionally, several other applications employing freeform surfaces have been demonstrated. Alvarez^{36} applied two complementary freeform optical components represented by third-order XY polynomials to realize a varying-focal-length lens system. Smilie et al.^{37} fabricated a germanium Alvarez lens for infrared imaging under conditions of modern precision optical manufacturing. Hicks^{38} developed a ray-tracing method to obtain a discrete data point cloud for a freeform reflector, which was used to design a 45-deg field of view and a lower-distortion driver-side mirror to effectively increase the driver’s field of view. Zhu et al.^{39} presented a direct design method to obtain a discrete data point cloud, which was fitted by XY polynomials for further optimization of an $F$-theta freeform scanning lens. Its field of view was $\sim 120\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$, and the scanning error was lower than $1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$.

The current applications of optical freeform surfaces show that there are mainly two strategies for characterizing the freeform surface. First, the freeform surface is expressed as explicit mathematical functions in the freeform optical design by directly optimizing the coefficients. Most of the commonly used freeform functions have been integrated into commercial optical design software. Second, a freeform surface is derived from discrete data points by specific design algorithms. Additionally, in freeform surface testing, the obtained freeform shape data are also included in a discrete data point cloud. They must be further fitted accurately by analytical or numerical functions for the next optimization or the final freeform surface estimation. Among the different function types for representing freeform surfaces, there is a particular impression about the XY-type polynomial freeform surface, which is frequently used in different applications. This can be illustrated briefly using the optical aberration theory. When the symmetry of an optical system is not considered, its aberration function can be approximated by expanding it as a combination of XY polynomials in a Taylor series for a specific field of view. In other words, each term of the XY polynomials corresponds to a specific aberration and can correct it. Similarly, the aberration function can be decomposed as a linear combination of Zernike polynomials, which is why Zernike polynomials are increasingly used. They provide more degrees of freedom and aberration correction capabilities in the freeform optical design. The advantages of freeform surfaces make them attractive to optical engineers for practical applications. In the following section, we describe the different freeform surface representation techniques in detail.

## 3.

## Mathematical Descriptions of Optical Freeform Surface

According to the different applications of freeform surfaces mentioned in Sec. 2, the common analytical functions used to represent optical freeform surfaces can be generally categorized as orthogonal polynomials and nonorthogonal functions. On the other hand, representation techniques employed to obtain an optical freeform surface from discrete data points are also important in direct freeform optics design and in freeform optical testing. For a freeform surface with large slope variation in local areas or large global gradient, hybrid or combined representation methods should be considered for characterizing the fine local features of a freeform surface. In this section, these representation techniques for freeform surfaces are discussed and analyzed.

## 3.1.

### Analytical Functions

## 3.1.1.

#### Orthogonal polynomials

Orthogonal polynomials,^{40} such as Zernike polynomials,^{41} are widely used to analyze optical surface deviations and wavefront aberrations owing to their elegant mathematical performance. Here, we describe mathematical orthogonality and different analytical orthogonal polynomials.

### Definition and properties of orthogonal polynomial

Polynomial functions ${P}_{m}(\rho )$, ${P}_{n}(\rho )$, and the weighted function $w(\rho )$ that are defined over a range $\rho \in [a,b]$ satisfy a relationship

then $P(\rho )$ is a weighted orthogonal polynomial, where ${\delta}_{mn}$ is the Kronecker delta function. The subscripts $m$ and $n$ are the nonnegative integers. When $m=n$, ${\delta}_{mn}=1$, and when $m\ne n$, ${\delta}_{mn}=0$. ${c}_{n}$ is the corresponding coefficient, which can be calculated by Eq. (2). When ${c}_{n}=1$, the polynomial is the orthonormal polynomialThe Zernike circle polynomial is a typical orthogonal and complete polynomial used to describe the aberration function, fit wavefront data, or represent a freeform surface. When the freeform surface is expanded as a linear combination of Zernike circle polynomials, Zernike circle polynomial has the following properties.

1. Because of the orthogonality and completeness of Zernike polynomials, the coefficient of one term is independent of the number of used terms.

2. In the expansion, the first term is the piston term, which indicates the mean value of the surface or wavefront function.

3. Except the first term, the mean value of the other terms is zero.

4. Except the first term, the variance of the surface or wavefront function is the sum of the squared coefficients of the other terms.

5. Except the first term, the standard deviation of each of the other terms is the corresponding coefficient value.

Obviously, other different-type orthogonal polynomials have similar properties to those of Zernike polynomials for representing an optical freeform surface.

### Two typical circular orthogonal polynomials

### Zernike circle orthogonal polynomials

The ordering number of the Zernike circle polynomial follows Noll’s ordering.^{42} A Zernike polynomial is expressed as

## (3)

$${Z}_{j}(\rho ,\varphi )={Z}_{n}^{m}(\rho ,\varphi )=\{\begin{array}{ll}\sqrt{n+1}{R}_{n}^{0}(\rho ),& m=0\\ \sqrt{2(n+1)}{R}_{n}^{m}(\rho )\mathrm{cos}\text{\hspace{0.17em}}m\varphi ,& m\ne 0\text{\hspace{0.17em}}(\text{even})\\ \sqrt{2(n+1)}{R}_{n}^{m}(\rho )\mathrm{sin}\text{\hspace{0.17em}}m\varphi ,& m\ne 0\text{\hspace{0.17em}}(\text{odd})\end{array},$$## (4)

$${R}_{n}^{m}(\rho )=\sum _{k=0}^{(n-m)/2}\frac{{(-1)}^{k}(n-k)!}{k!(\frac{n+m}{2}-k)!(\frac{n-m}{2}-k)!}{\rho}^{n-2k}.$$The higher-order Zernike polynomials can be obtained from the recurrence relations that they satisfy. The first 28 Zernike polynomials are listed in Table 1.

## Table 1

The first 28 Zernike circle polynomials Zj(ρ,ϕ)

n | m | j | Zj(ρ,ϕ) |
---|---|---|---|

0 | 0 | 1 | 1 |

1 | 1 | 2 | $2\rho \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\varphi $ |

1 | 1 | 3 | $2\rho \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\varphi $ |

2 | 0 | 4 | $\sqrt{3}(2{\rho}^{2}-1)$ |

2 | 2 | 5 | $\sqrt{6}{\rho}^{2}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\varphi $ |

2 | 2 | 6 | $\sqrt{6}{\rho}^{2}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\varphi $ |

3 | 1 | 7 | $\sqrt{8}(3{\rho}^{3}-2\rho )\mathrm{sin}\text{\hspace{0.17em}}\varphi $ |

3 | 1 | 8 | $\sqrt{8}(3{\rho}^{3}-2\rho )\mathrm{cos}\text{\hspace{0.17em}}\varphi $ |

3 | 3 | 9 | $\sqrt{8}{\rho}^{3}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}3\varphi $ |

3 | 3 | 10 | $\sqrt{8}{\rho}^{3}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}3\varphi $ |

4 | 0 | 11 | $\sqrt{5}(6{\rho}^{4}-6{\rho}^{2}+1)$ |

4 | 2 | 12 | $\sqrt{10}(4{\rho}^{4}-3{\rho}^{2})\mathrm{cos}\text{\hspace{0.17em}}2\varphi $ |

4 | 2 | 13 | $\sqrt{10}(4{\rho}^{4}-3{\rho}^{2})\mathrm{sin}\text{\hspace{0.17em}}2\varphi $ |

4 | 4 | 14 | $\sqrt{10}{\rho}^{4}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}4\varphi $ |

4 | 4 | 15 | $\sqrt{10}{\rho}^{4}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}4\varphi $ |

5 | 1 | 16 | $\sqrt{12}(10{\rho}^{5}-12{\rho}^{3}+3\rho )\mathrm{cos}\text{\hspace{0.17em}}\varphi $ |

5 | 1 | 17 | $\sqrt{12}(10{\rho}^{5}-12{\rho}^{3}+3\rho )\mathrm{sin}\text{\hspace{0.17em}}\varphi $ |

5 | 3 | 18 | $\sqrt{12}(5{\rho}^{5}-4{\rho}^{3})\mathrm{cos}\text{\hspace{0.17em}}3\varphi $ |

5 | 3 | 19 | $\sqrt{12}(5{\rho}^{5}-4{\rho}^{3})\mathrm{sin}\text{\hspace{0.17em}}3\varphi $ |

5 | 5 | 20 | $\sqrt{12}{\rho}^{5}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}5\varphi $ |

5 | 5 | 21 | $\sqrt{12}{\rho}^{5}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}5\varphi $ |

6 | 0 | 22 | $\sqrt{7}(20{\rho}^{6}-30{\rho}^{4}+12{\rho}^{2}-1)$ |

6 | 2 | 23 | $\sqrt{14}(15{\rho}^{6}-20{\rho}^{4}+6{\rho}^{2})\mathrm{sin}\text{\hspace{0.17em}}2\varphi $ |

6 | 2 | 24 | $\sqrt{14}(15{\rho}^{6}-20{\rho}^{4}+6{\rho}^{2})\mathrm{cos}\text{\hspace{0.17em}}2\varphi $ |

6 | 4 | 25 | $\sqrt{14}(6{\rho}^{6}-5{\rho}^{4})\mathrm{sin}\text{\hspace{0.17em}}4\varphi $ |

6 | 4 | 26 | $\sqrt{14}(6{\rho}^{6}-5{\rho}^{4})\mathrm{cos}\text{\hspace{0.17em}}4\varphi $ |

6 | 6 | 27 | $\sqrt{14}{\rho}^{6}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}6\varphi $ |

6 | 6 | 28 | $\sqrt{14}{\rho}^{6}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}6\varphi $ |

As discussed in Sec. 2, the Zernike polynomial is widely used. Fuerschbach et al.^{9} designed a compact off-axis three-mirror long-wavelength infrared imaging system using fringe Zernike polynomials to further improve its imaging performance. In theory, any surface or wavefront, regardless of how complex it is, can be accurately represented by applying a sufficient number of Zernike polynomial terms. However, in practice, only a finite number of terms are usually used for freeform surface analysis.

### Q-type orthogonal polynomials

To overcome the numerical deficiency of even-order polynomials and improve the manufacturability and testability of aspheres, Forbes proposed ${Q}_{n}^{\mathrm{con}}({u}^{2})$ polynomials for representing strong surfaces with large aspheric deviations and ${Q}_{n}^{bfs}({u}^{2})$ polynomials for characterizing mild surfaces with constrained slopes based on Jacobi polynomials.^{43}44.45.^{–}^{46} Furthermore, a ${Q}_{n}^{m}({u}^{2})$ polynomial was proposed for characterizing optical freeform surfaces.^{47} The ${Q}_{n}^{m}({u}^{2})$ polynomial has a similar construction to the Zernike polynomial and represents the deviation between the freeform surface and the best-fit sphere along the normal direction. The sag of the freeform surface along the $z$-axis is $z=f(\rho ,\varphi )$

## (5)

$$z=\frac{{c}_{bfs}{\rho}^{2}}{1+\sqrt{1-{c}_{bfs}^{2}{\rho}^{2}}}+\frac{1}{\sqrt{1-{c}_{bfs}^{2}{\rho}^{2}}}[{u}^{2}(1-{u}^{2})\sum _{n=0}^{N}{a}_{n}^{0}{Q}_{n}^{0}({u}^{2})\phantom{\rule{0ex}{0ex}}+\sum _{m=1}^{M}{u}^{m}\sum _{n=0}^{N}({a}_{n}^{m}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}m\varphi +{b}_{n}^{m}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}m\varphi ){Q}_{n}^{m}({u}^{2})].$$## (6)

$${J}_{n}^{m}(t)=\{\begin{array}{ll}1-\frac{t}{2},& m=1\text{\hspace{0.17em}\hspace{0.17em}}\text{and}\text{\hspace{0.17em}\hspace{0.17em}}n=1\\ \frac{{(-1)}^{n}(2n)!!}{2(2n-1)!!}{J}_{n}^{(-\frac{3}{2},m-\frac{3}{2})}(2t-1)& \text{otherwise}\end{array}.$$## Table 2

The first four-terms of Qn6(t) polynomials.

n | Qn6(t) |
---|---|

0 | $\frac{8}{9\sqrt{7}}$ |

1 | $-\frac{8(-77+72t)}{9\sqrt{1397}}$ |

2 | $\frac{8(73645-159432t+85344{t}^{2})}{3\sqrt{542196655}}$ |

3 | $-\frac{8(-1858285+6903936t-8322912{t}^{2}+3275520{t}^{3})}{9\sqrt{4280409445}}$ |

To obtain the higher-order Q-type polynomials, a recurrence relation can be applied to improve the computation efficiency.

### Different orthogonal polynomials with differently shaped apertures

For the circular aperture or pupil of the optical system, Zernike polynomials provide a superior performance. In modern optical engineering, more optical systems with noncircular pupils or optical components with noncircular boundaries are being developed. For example, rectangular or square optical elements are applied in anamorphic or high-power laser optical systems.^{48} Many optical elements with square apertures are used in the U.S. National Ignition Facility and the China Inertial Confinement Fusion laser driver. In large telescopes, the giant reflector is stitched by hundreds of hexagonal mirrors. Freeform optical elements with noncircular apertures are also implemented in practical applications.^{49}

For freeform surface analysis in freeform optical elements with differently shaped apertures, the corresponding orthogonal polynomials must be applied to the similarly shaped apertures. Owing to the orthogonality and completeness of Zernike polynomials, a freeform surface with noncircular aperture can also be expressed as a linear combination of the corresponding orthogonal polynomials over the same aperture. Using Zernike polynomials as the basis functions, orthogonal polynomials (such as Zernike square polynomials, Zernike rectangular polynomials, Zernike hexagonal polynomials, Zernike elliptical polynomials, and Zernike annular polynomials) over the corresponding apertures can be derived by the Gram–Schmidt orthogonalization method.^{50} Tatian,^{51} Barakat,^{52} Mahajan and Dai,^{50} Swantner and Chow,^{53} Hasan and Shaker,^{54} Díaz and Mahajan,^{55} and Díaz and Navarro^{56} reported analytical Zernike orthogonal polynomials with differently shaped apertures. Ferreira et al.^{57} proposed a rigorous and powerful theoretical framework to obtain the orthogonal basis with a conicoid first mode for surface specification. Liu et al.^{58} applied two-dimensional (2-D) Chebyshev polynomials to characterize the “W”-shaped freeform optical elements. Mahajan^{59} analyzed the aberration of an anamorphic optical system using 2-D Legendre polynomials. These two orthogonal polynomials were used for square or rectangular apertures. The differences between 2-D Chebyshev and Legendre polynomials for representing complex freeform surfaces were compared comprehensively by Ye et al.,^{60} who elaborated on the mathematical expressions of several typical orthogonal polynomials for representing freeform surfaces.

### Gram–Schmidt orthogonalization method deriving noncircular orthogonal polynomials

For freeform surfaces with noncircular apertures, such as square, rectangular, elliptical, annular, sector, and hexagonal apertures, the corresponding orthogonal polynomials can be derived by the Gram–Schmidt orthogonalization method,^{50} which employs a set of complete orthogonal polynomials as basis functions and uses an iterative transformation method.

For simplicity, a Zernike circle polynomial ${Z}_{j}(x,y)$ in the Cartesian coordinate system is written as ${Z}_{j}$. The subscript $j$ is a positive integer starting from 1 and the ordering number, ${F}_{j}$ denotes the corresponding noncircular orthogonal polynomial, and ${G}_{j}$ is the intermediate transformation polynomial. Generally, the first term is

Furthermore, if $G$, $Z$, and $F$ satisfy the recurrence relation ${G}_{j+1}={Z}_{j+1}+\sum _{k=1}^{j}{c}_{j+1,k}{F}_{k}$, then the corresponding noncircular normalized orthogonal polynomial $F$ is

## (8)

$${F}_{j+1}=\frac{{G}_{j+1}}{\Vert {G}_{j+1}\Vert}=\frac{{G}_{j+1}}{{\left(\frac{1}{A}{\int}_{\text{aperture}}{G}_{j+1}^{2}\mathrm{d}x\text{\hspace{0.17em}}\mathrm{d}y\right)}^{1/2}}\mathrm{.}$$The transformation coefficient ${c}_{j+1,k}$ is

## (9)

$${c}_{j+1,k}=-\frac{1}{A}{\int}_{\text{aperture}}{Z}_{j+1}{F}_{k}\mathrm{d}x\text{\hspace{0.17em}}\mathrm{d}y,$$Using recurrence iteration based on the Gram–Schmidt orthogonalization method, analytical orthogonal polynomials for differently shaped apertures can be derived. These polynomials can be used for freeform surface representation, wavefront estimation, and aberration analysis.

### Two typical Jacobi square orthogonal polynomials

In the orthogonal polynomial used for the square aperture derived by the Gram–Schmidt orthogonalization method, the $x$ and $y$ variables of each term are mixed together. The 2-D Chebyshev and Legendre polynomials are two square orthogonal polynomials with $x$ and $y$ variables, respectively, which can be applied to characterize anamorphic freeform surfaces.

### Two-dimensional Chebyshev polynomials

The 2-D Chebyshev polynomials are defined by the products of one-dimensional first-kind Chebyshev polynomials in the $x$- and $y$-dimensions. The Chebyshev polynomials ${T}_{n}(x)$ of the first kind^{61} in the $x$-dimension are orthogonal in the range $[-\mathrm{1,1}]$ with a weighted function $w(x)=1/\sqrt{1-{x}^{2}}$, as shown in

## (10)

$${\int}_{-1}^{1}\frac{{T}_{n}(x){T}_{m}(x)}{\sqrt{1-{x}^{2}}}\mathrm{d}x=\{\begin{array}{ll}\frac{1}{2}\pi {\delta}_{nm}& \text{for}\text{\hspace{0.17em}\hspace{0.17em}}n\ne 0,m\ne 0\\ \pi & \text{for}\text{\hspace{0.17em}\hspace{0.17em}}n=m=0\end{array},$$## (12)

$${\int}_{-1}^{1}{\int}_{-1}^{1}{C}_{j}(x,y){C}_{{j}^{\prime}}(x,y)\frac{\mathrm{d}x\text{\hspace{0.17em}}\mathrm{d}y}{\sqrt{1-{x}^{2}}\sqrt{1-{y}^{2}}}=\{\begin{array}{cc}0& \text{for}\text{\hspace{0.17em}\hspace{0.17em}}j\ne {j}^{\prime}\\ K& \text{for}\text{\hspace{0.17em}\hspace{0.17em}}j={j}^{\prime}\end{array},$$## (13)

$$K=\{\begin{array}{ll}{\pi}^{2}& \text{for}\text{\hspace{0.17em}\hspace{0.17em}}n=m=0\\ {\pi}^{2}/4& \text{for}\text{\hspace{0.17em}\hspace{0.17em}}n=m\ne 0\\ {\pi}^{2}/2& \text{for}\text{\hspace{0.17em}\hspace{0.17em}}\text{else}\end{array}.$$The normalized 2-D Chebyshev polynomial ${C}_{Nj}$ can be derived by Eq. (14). The first 28 2-D Chebyshev polynomials are shown in Table 3.

## (14)

$${C}_{Nj}=\frac{{C}_{j}/\sqrt{K}}{\Vert {C}_{j}/\sqrt{K}\Vert}=\frac{{C}_{j}/\sqrt{K}}{{\left[\frac{1}{4}{\int}_{-1}^{1}{\int}_{-1}^{1}{({C}_{j}/\sqrt{K})}^{2}\mathrm{d}x\text{\hspace{0.17em}}\mathrm{d}y\right]}^{1/2}}.$$## Table 3

The first 28 2-D Chebyshev polynomials Cj(x,y).

Polynomial order | Cj(x,y) |
---|---|

0 | ${C}_{1}=1$ |

1 | ${C}_{2}=x$ |

1 | ${C}_{3}=y$ |

2 | ${C}_{4}=2{x}^{2}-1$ |

2 | ${C}_{5}=xy$ |

2 | ${C}_{6}=2{y}^{2}-1$ |

3 | ${C}_{7}=4{x}^{3}-3x$ |

3 | ${C}_{8}=(2{x}^{2}-1)y$ |

3 | ${C}_{9}=x(2{y}^{2}-1)$ |

3 | ${C}_{10}=4{y}^{3}-3y$ |

4 | ${C}_{11}=8{x}^{4}-8{x}^{2}+1$ |

4 | ${C}_{12}=(4{x}^{3}-3x)y$ |

4 | ${C}_{13}=(2{x}^{2}-1)(2{y}^{2}-1)$ |

4 | ${C}_{14}=x(4{y}^{3}-3y)$ |

4 | ${C}_{15}=8{y}^{4}-8{y}^{2}+1$ |

5 | ${C}_{16}=16{x}^{5}-20{x}^{3}+5x$ |

5 | ${C}_{17}=(8{x}^{4}-8{x}^{2}+1)y$ |

5 | ${C}_{18}=(4{x}^{3}-3x)(2{y}^{2}-1)$ |

5 | ${C}_{19}=(2{x}^{2}-1)(4{y}^{3}-3y)$ |

5 | ${C}_{20}=x(8{y}^{4}-8{y}^{2}+1)$ |

5 | ${C}_{21}=16{y}^{5}-20{y}^{3}+5y$ |

6 | ${C}_{22}=32{x}^{6}-48{x}^{4}+18{x}^{2}-1$ |

6 | ${C}_{23}=(16{x}^{5}-20{x}^{3}+5x)y$ |

6 | ${C}_{24}=(8{x}^{4}-8{x}^{2}+1)(2{y}^{2}-1)$ |

6 | ${C}_{25}=(4{x}^{3}-3x)(4{y}^{3}-3y)$ |

6 | ${C}_{26}=(2{x}^{2}-1)(8{y}^{4}-8{y}^{2}+1)$ |

6 | ${C}_{27}=x(16{y}^{5}-20{y}^{3}+5y)$ |

6 | ${C}_{28}=32{y}^{6}-48{y}^{4}+18{y}^{2}-1$ |

### Two-dimensional Legendre polynomials

The 2-D Legendre polynomials are obtained in a similar way to 2-D Chebyshev polynomials. The Legendre polynomials^{62} ${P}_{n}(x)$ in the $x$-dimension are orthogonal over the interval $[-\mathrm{1,1}]$ with a weighted function $w(x)=1$

To obtain ${P}_{n}(y)$, we substitute the $x$ variable in ${P}_{n}(x)$ with the $y$ variable. Thus, the 2-D Legendre polynomials ${L}_{j}(x,y)$ can be obtained from the product of ${P}_{n}(x)$ and ${P}_{m}(y)$

The orthogonality of 2-D Legendre polynomials is described by

## (17)

$$\frac{1}{4}{\int}_{-1}^{1}{\int}_{-1}^{1}{L}_{j}(x,y){L}_{{j}^{\prime}}(x,y)\mathrm{d}x\text{\hspace{0.17em}}\mathrm{d}y={\delta}_{j{j}^{\prime}}.$$The first 28 2-D Legendre polynomials are listed in Table 4.

## Table 4

The first 28 2-D Legendre polynomials Lj(x,y).

Polynomial order | Lj(x,y) |
---|---|

0 | ${L}_{1}=1$ |

1 | ${L}_{2}=\sqrt{3}x$ |

1 | ${L}_{3}=\sqrt{3}y$ |

2 | ${L}_{4}=\sqrt{5}(3{x}^{2}-1)/2$ |

2 | ${L}_{5}=3xy$ |

2 | ${L}_{6}=\sqrt{5}(3{y}^{2}-1)/2$ |

3 | ${L}_{7}=\sqrt{7}(5{x}^{3}-3x)/2$ |

3 | ${L}_{8}=\sqrt{15}(3{x}^{2}-1)y/2$ |

3 | ${L}_{9}=\sqrt{15}x(3{y}^{2}-1)/2$ |

3 | ${L}_{10}=\sqrt{7}(5{y}^{3}-3y)/2$ |

4 | ${L}_{11}=3(35{x}^{4}-30{x}^{2}+3)/8$ |

4 | ${L}_{12}=\sqrt{21}(5{x}^{3}-3x)y/2$ |

4 | ${L}_{13}=5(3{x}^{2}-1)(3{y}^{2}-1)/4$ |

4 | ${L}_{14}=\sqrt{21}x(5{y}^{3}-3y)/2$ |

4 | ${L}_{15}=3(35{y}^{4}-30{y}^{2}+3)/8$ |

5 | ${L}_{16}=\sqrt{11}(63{x}^{5}-70{x}^{3}+15x)/8$ |

5 | ${L}_{17}=3\sqrt{3}(35{x}^{4}-30{x}^{2}+3)y/8$ |

5 | ${L}_{18}=\sqrt{35}(5{x}^{3}-3x)(3{y}^{2}-1)/4$ |

5 | ${L}_{19}=\sqrt{35}(3{x}^{2}-1)(5{y}^{3}-3y)/4$ |

5 | ${L}_{20}=3\sqrt{3}x(35{y}^{4}-30{y}^{2}+3)/8$ |

5 | ${L}_{21}=\sqrt{11}(63{y}^{5}-70{y}^{3}+15y)/8$ |

6 | ${L}_{22}=\sqrt{13}(231{x}^{6}-315{x}^{4}+105{x}^{2}-5)/16$ |

6 | ${L}_{23}=\sqrt{33}(63{x}^{5}-70{x}^{3}+15x)y/8$ |

6 | ${L}_{24}=3\sqrt{5}(35{x}^{4}-30{x}^{2}+3)(3{y}^{2}-1)/16$ |

6 | ${L}_{25}=7(5{x}^{3}-3x)(5{y}^{3}-3y)/4$ |

6 | ${L}_{26}=3\sqrt{5}(3{x}^{2}-1)(35{y}^{4}-30{y}^{2}+3)/16$ |

6 | ${L}_{27}=\sqrt{33}x(63{y}^{5}-70{y}^{3}+15y)/8$ |

6 | ${L}_{28}=\sqrt{13}(231{y}^{6}-315{y}^{4}+105{y}^{2}-5)/16$ |

## 3.1.2.

#### Nonorthogonal functions

In addition to the analytical orthogonal polynomials, there are several nonorthogonal functions applied to characterize freeform surfaces, such as XY-type polynomials, spline functions, and radial basis functions.

### XY-type polynomials

In XY-type polynomials, the freeform surface is represented as a base surface with a linear combination of multiple ${x}^{m}{y}^{n}$ monomial terms. The base surface is usually a conic or anamorphic surface, as shown in

## (18)

$$z(x,y)=\frac{{c}_{x}{x}^{2}+{c}_{y}{y}^{2}}{1+\sqrt{1-({k}_{x}+1){c}_{x}^{2}{x}^{2}-({k}_{y}+1){c}_{y}^{2}{y}^{2}}}+\sum _{j=2}^{J}{a}_{j}{x}^{m}{y}^{n},\phantom{\rule[-0.0ex]{2em}{0.0ex}}j=\frac{{(m+n)}^{2}+m+3n}{2}+1,$$This type of freeform surface is not available in optical design software. However, it can be defined in such software using the feature “user-defined surface” for further optimization of the corresponding freeform optical system.^{63} Additionally, the off-axis part of the XY-type polynomial is usually used to represent freeform surfaces.

### Spline surface

A spline surface is very flexible in computer-aided design and freeform surface construction. The NURBS surface is increasingly being applied to optical design. Its mathematical expression is

## (19)

$$z(u,v)=\frac{\sum _{i=0}^{n}\sum _{j=0}^{m}{w}_{i,j}{B}_{i,k}(u){B}_{j,l}(v){P}_{i,j}}{\sum _{i=0}^{n}\sum _{j=0}^{m}{w}_{i,j}{B}_{i,k}(u){B}_{j,l}(v)},$$## (20)

$$\{\begin{array}{ll}k=0,& {B}_{i,0}(u)=\{\begin{array}{ll}1,& \text{if}\text{\hspace{0.17em}\hspace{0.17em}}{u}_{i}\le u<{u}_{i+1}\text{\hspace{0.17em}\hspace{0.17em}}\text{and}\text{\hspace{0.17em}\hspace{0.17em}}{u}_{i}<{u}_{i+1}\\ 0,& \text{otherwise}\end{array}\\ k\ge 1,& {B}_{i,k}(u)=\frac{u-{u}_{i}}{{u}_{i+k}-{u}_{i}}{B}_{i,k-1}(u)+\frac{{u}_{i+k+1}-u}{{u}_{i+k+1}-{u}_{i+1}}{B}_{i+1,k-1}(u)\end{array},$$## (21)

$$\{\begin{array}{ll}l=0,& {B}_{j,0}(v)=\{\begin{array}{ll}1,& \text{if}\text{\hspace{0.17em}\hspace{0.17em}}{v}_{j}\le v<{v}_{j+1}\text{\hspace{0.17em}\hspace{0.17em}}\text{and}\text{\hspace{0.17em}\hspace{0.17em}}{v}_{j}<{v}_{j+1}\\ 0,& \text{otherwise}\end{array}\\ l\ge 1,& {B}_{j,l}(v)=\frac{v-{v}_{j}}{{v}_{j+l}-{v}_{j}}{B}_{j,l-1}(v)+\frac{{v}_{j+l+1}-v}{{v}_{j+l+1}-{v}_{j+1}}{B}_{j+1,l-1}(v)\end{array}.$$The local feature of a NURBS surface can be controlled by the control function and weighted factor. The freeform shape can be changed flexibly. Chrisp et al.^{64}^{,}^{65} designed an imaging freeform optical system using NURBS freeform surfaces. The imaging performance of the system was improved compared with that obtained when using Zernike polynomial surfaces. Radial and toroidal NURBS surfaces have been integrated into commercial optical design software, which gives optical designers more options. The next step should be improving the ray-tracing speed and efficiency when using NURBS functions for characterizing freeform surfaces.

### Radial basis function

A linear combination of radial basis functions with the conic surface forms another freeform surface type, as described by

## (22)

$$z(x,y)=\frac{c({x}^{2}+{y}^{2})}{1+\sqrt{1-(k+1){c}^{2}({x}^{2}+{y}^{2})}}+\sum _{n=1}^{N}{w}_{n}{\mathrm{\varphi}}_{n},$$Eq. (23) shows that the supported center ${\mathbf{p}}_{n}$ can be changed in the supported region. Thus, the radial basis function is multicentric. Cakmakci et al.^{27} designed a head-worn display using a Gaussian radial basis function freeform surface, which was integrated into optical design software by applying the “user-defined surface” feature. Moreover, the local feature of the freeform surface can be controlled by the radial basis function.

## 3.2.

### Representation Techniques for Deriving an Optical Freeform Surface from Discrete Data Points

## 3.2.1.

#### Numerical orthogonal polynomials

The analytical orthogonal polynomials with differently shaped apertures derived by the Gram–Schmidt orthogonalization method are, in theory, only orthogonal over their continuous domain of definition. Generally, the obtained data are the discrete data points, where the orthogonality of the analytical polynomials may be degraded. Thus, a set of orthogonal polynomials that can adapt discrete data points are necessary. On the other hand, for a freeform surface with a complex aperture, applying the Gram–Schmidt orthogonalization method with iterations is rather tedious. Malacara et al.^{66} presented discrete orthogonal polynomials for wavefront analysis with a circular aperture. Dai and Mahajan^{67} proposed a noniterative and fast algorithm to obtain analytical orthogonal polynomials with arbitrarily shaped aperture. Furthermore, Ye et al.^{49} presented numerical orthogonal polynomials based on the matrix transformation method for fitting a discrete data point cloud of a freeform surface.

Initially, a numerical orthogonal polynomial is expressed as a linear combination of Zernike polynomials owing to its orthogonality and completeness, as shown in

where ${F}_{l}({x}_{n},{y}_{n})$ is the numerical orthogonal polynomial, the subscript $l$ is the ordering number, and $({x}_{n},{y}_{n})$ is the point coordinate in the effective region. The number of data points in the effective region is $N$, where $n=\mathrm{1,2},\dots ,N$. ${M}_{lj}$ is the transformation coefficient, and $J$ is the number of terms of the Zernike polynomial ${Z}_{j}({x}_{n},{y}_{n})$.For a data point, Eq. (24) can be expanded as

## (25)

$$\{\begin{array}{l}{F}_{1}({x}_{n},{y}_{n})={M}_{11}{Z}_{1}({x}_{n},{y}_{n})+{M}_{12}{Z}_{2}({x}_{n},{y}_{n})+\cdots +{M}_{1j}{Z}_{j}({x}_{n},{y}_{n})+\cdots +{M}_{1J}{Z}_{J}({x}_{n},{y}_{n})\\ {F}_{2}({x}_{n},{y}_{n})={M}_{21}{Z}_{1}({x}_{n},{y}_{n})+{M}_{22}{Z}_{2}({x}_{n},{y}_{n})+\cdots +{M}_{2j}{Z}_{j}({x}_{n},{y}_{n})+\cdots +{M}_{2J}{Z}_{J}({x}_{n},{y}_{n})\\ \vdots \\ {F}_{j}({x}_{n},{y}_{n})={M}_{j1}{Z}_{1}({x}_{n},{y}_{n})+{M}_{j2}{Z}_{2}({x}_{n},{y}_{n})+\cdots +{M}_{jj}{Z}_{j}({x}_{n},{y}_{n})+\cdots +{M}_{jJ}{Z}_{J}({x}_{n},{y}_{n})\\ \vdots \\ {F}_{J}({x}_{n},{y}_{n})={M}_{J1}{Z}_{1}({x}_{n},{y}_{n})+{M}_{J2}{Z}_{2}({x}_{n},{y}_{n})+\cdots +{M}_{Jj}{Z}_{j}({x}_{n},{y}_{n})+\cdots +{M}_{JJ}{Z}_{J}({x}_{n},{y}_{n})\end{array}.$$The transform matrix $\mathbf{M}$ is

## (26)

$$\mathbf{M}=\left[\begin{array}{cccccc}{M}_{11}& {M}_{12}& \cdots & {M}_{1j}& \cdots & {M}_{1J}\\ {M}_{21}& {M}_{22}& \cdots & {M}_{2j}& \cdots & {M}_{2J}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {M}_{j1}& {M}_{j2}& \cdots & {M}_{jj}& \cdots & {M}_{jJ}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {M}_{J1}& {M}_{J2}& \cdots & {M}_{Jj}& \cdots & {M}_{JJ}\end{array}\right].$$Therefore, Eq. (25) can be described in the matrix form

## (27)

$$[\begin{array}{cccccc}{F}_{1}({x}_{n},{y}_{n})& {F}_{2}({x}_{n},{y}_{n})& \cdots & {F}_{j}({x}_{n},{y}_{n})& \cdots & {F}_{J}({x}_{n},{y}_{n})\end{array}]\phantom{\rule{0ex}{0ex}}=[\begin{array}{cccccc}{Z}_{1}({x}_{n},{y}_{n})& {Z}_{2}({x}_{n},{y}_{n})& \cdots & {Z}_{j}({x}_{n},{y}_{n})& \cdots & {Z}_{J}({x}_{n},{y}_{n})\end{array}]\left[\begin{array}{cccccc}{M}_{11}& {M}_{21}& \cdots & {M}_{j1}& \cdots & {M}_{J1}\\ {M}_{12}& {M}_{22}& \cdots & {M}_{j2}& \cdots & {M}_{J2}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {M}_{1j}& {M}_{2j}& \cdots & {M}_{jj}& \cdots & {M}_{Jj}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {M}_{1J}& {M}_{2J}& \cdots & {M}_{jJ}& \cdots & {M}_{JJ}\end{array}\right].$$For all the data points in the effective region, it can be

## (28)

$$\left[\begin{array}{cccccc}{F}_{1}({x}_{1},{y}_{1})& {F}_{2}({x}_{1},{y}_{1})& \cdots & {F}_{j}({x}_{1},{y}_{1})& \cdots & {F}_{J}({x}_{1},{y}_{1})\\ {F}_{1}({x}_{2},{y}_{2})& {F}_{2}({x}_{2},{y}_{2})& \cdots & {F}_{j}({x}_{2},{y}_{2})& \cdots & {F}_{J}({x}_{2},{y}_{2})\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {F}_{1}({x}_{n},{y}_{n})& {F}_{2}({x}_{n},{y}_{n})& \cdots & {F}_{j}({x}_{n},{y}_{n})& \cdots & {F}_{J}({x}_{n},{y}_{n})\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {F}_{1}({x}_{N},{y}_{N})& {F}_{2}({x}_{N},{y}_{N})& \cdots & {F}_{j}({x}_{N},{y}_{N})& \cdots & {F}_{J}({x}_{N},{y}_{N})\end{array}\right]\phantom{\rule{0ex}{0ex}}=\left[\begin{array}{cccccc}{Z}_{1}({x}_{1},{y}_{1})& {Z}_{2}({x}_{1},{y}_{1})& \cdots & {Z}_{j}({x}_{1},{y}_{1})& \cdots & {Z}_{J}({x}_{1},{y}_{1})\\ {Z}_{1}({x}_{2},{y}_{2})& {Z}_{2}({x}_{2},{y}_{2})& \cdots & {Z}_{j}({x}_{2},{y}_{2})& \cdots & {Z}_{J}({x}_{2},{y}_{2})\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {Z}_{1}({x}_{n},{y}_{n})& {Z}_{2}({x}_{n},{y}_{n})& \cdots & {Z}_{j}({x}_{n},{y}_{n})& \cdots & {Z}_{J}({x}_{n},{y}_{n})\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {Z}_{1}({x}_{N},{y}_{N})& {Z}_{2}({x}_{N},{y}_{N})& \cdots & {Z}_{j}({x}_{N},{y}_{N})& \cdots & {Z}_{J}({x}_{N},{y}_{N})\end{array}\right]\left[\begin{array}{cccccc}{M}_{11}& {M}_{21}& \cdots & {M}_{j1}& \cdots & {M}_{J1}\\ {M}_{12}& {M}_{22}& \cdots & {M}_{j2}& \cdots & {M}_{J2}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {M}_{1j}& {M}_{2j}& \cdots & {M}_{jj}& \cdots & {M}_{Jj}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {M}_{1\text{\hspace{0.17em}\hspace{0.17em}}J}& {M}_{2\text{\hspace{0.17em}\hspace{0.17em}}J}& \cdots & {M}_{jJ}& \cdots & {M}_{JJ}\end{array}\right].$$The derivation process of a numerical matrix $\mathbf{F}$ formed by the numerical orthogonal polynomial ${F}_{l}({x}_{n},{y}_{n})$ is in general. According to its practical use, a numerical orthogonal polynomial can also be decomposed as a linear combination of other orthogonal polynomials with completeness. The numerical orthogonal polynomial can be applied to represent the discrete data point cloud of a freeform surface as well as a freeform surface with a general shaped aperture.

## 3.2.2.

#### Other representation techniques for handling discrete data points

In the freeform optical design, the lack of an original reference freeform optical configuration is relatively common. The optimization and design methods have been gradually developed. The direct design method for freeform surface is remarkable. It can be categorized as the aplanatic design method,^{68} partial differential equation method,^{69} SMS method,^{70} and iterative construction method.^{71} The direct design method obtains discrete data points for the freeform surface, and the discrete data point cloud must be fitted accurately by analytical or numerical functions for further optimization. In freeform optical testing, the obtained freeform shape data also form a discrete data point cloud. In theory, when the sampling dataset is sufficiently large, analytical orthogonal polynomials with the corresponding shape of aperture can be used approximately. Additionally, for a freeform surface with a complex aperture shape, the corresponding analytical orthogonal polynomial is usually difficult to obtain. Malacara et al.,^{66} Dai and Mahajan,^{67} Ye et al.,^{49} and Hilbig et al.^{72} have presented different methods to overcome this issue.

On the other hand, freeform surface or wavefront estimation from a slope measurement is an alternative approach used in modern optical testing, such as Shack–Hartmann sensing and phase measuring deflectometry. The measured data are the gradient-related discrete data point cloud. The methods used to derive a freeform surface or wavefront from its slope can be classified as zonal^{73} and modal^{74} methods. Using the zonal method, Southwell^{73} analyzed the differences of different data sampling types. Li et al.^{75} applied higher-order truncation errors based on the Taylor series expansion to increase the estimation accuracy in the Southwell data sampling type. Huang and Asundi^{76} used an iterative compensation method. For freeform surfaces with complex aperture shapes, the sag data of the freeform surface derived from the slope could be treated with the discrete Fourier transform^{77} or discrete cosine transform-based methods.^{78} Zou and Rolland^{79} presented an iterative zonal estimation method with the Gerchberg iteration^{80} to obtain a freeform surface for general shaped pupils. However, the zonal method is limited by the shape of the aperture, and the iterative method is not very efficient for dynamically varying aperture shapes.

In the modal method, the freeform surface or wavefront is decomposed as a linear combination of basis functions. The final coefficients characterizing the freeform surface are obtained from the measured slope data. Cubalchini^{74} proposed using the gradient function of Zernike polynomials as the basis function for fitting the measured gradient data. The drawback of this approach was that the gradient function of a Zernike polynomial is not an orthogonal function, which caused instability in the coefficients. Thus, a set of polynomials whose gradient functions are orthogonal over the measured region must be selected. Zhao and Burge^{81} proposed orthonormal vector polynomials in a unit circle. Obviously, for noncircular apertures, Zhao’s method could be limited. Mochi and Goldberg^{82} reported an iterative method for obtaining an orthogonal basis function to reconstruct a wavefront from its gradient. In this method, a finite set of 2-D polynomials that are expected to describe the measured wavefront over the testing domain must be initially selected. However, obtaining an orthogonal function over the aperture is not very easy, especially for complex-shaped apertures. Ye et al.^{83} presented a numerical orthogonal transformation method to obtain numerical orthogonal gradient polynomials for directly representing the measured gradient data of a freeform surface with a general shaped aperture. This method was noniterative and efficient and could also be applied to the slope-based freeform surface or wavefront sensing with the dynamically varying aperture shapes.

## 3.3.

### Representation Methods for Freeform Surfaces with Strong Slope Variation

Regarding freeform surfaces with strong slope variations, characterizing fine local features is a key issue for freeform surface estimation. When using only finite polynomial terms, the representation accuracy could be limited. Kaya et al.^{84} performed a comparative assessment between Zernike polynomials and Q-type polynomials for precisely characterizing a freeform surface. To fit an asymmetric local feature with sub-nanometer accuracy, the number of polynomial terms would reach thousands. The calculation of higher-order Zernike or Q-type polynomials is relatively tedious. Trevino et al.^{85} selected the first-class Bessel circular functions for characterizing a complex corneal surface. Svechnikov et al.^{86} discussed the lateral resolving capacity of circular Zernike polynomials and the accuracy of representing complex freeform surfaces, such as Gaussian and Gaussian-like surfaces. When more polynomial terms are applied, the efficiency of the freeform surface estimation is noticeably reduced.

Therefore, the hybrid or combination reconstruction method could be an alternative approach to freeform surface fitting with finite polynomials and high accuracy. The concept of the hybrid representation method is similar to that of subaperture stitching testing for large-aperture optical surfaces. The entire aperture of a freeform surface with strong slope variation is decomposed into multiple overlapping subapertures. Every local surface corresponding to the subaperture can be reconstructed by Zernike polynomials or other orthogonal polynomials. Then, the entire freeform surface is synthetized by all the local surfaces using specific algorithms. Espinosa et al.^{87} investigated a combination of zonal and modal fitting methods for evaluating an irregular corneal surface. Its zonal fitting by Zernike circle polynomial was implemented over a square subaperture rather than over a circular subaperture. Kaya and Rolland^{88} proposed a hybrid method combining local fitting in each subaperture and using a radial basis function as a global approximant for obtaining the entire surface. Ye et al.^{89} presented a combination method for representing complex freeform surfaces. Every local surface was reconstructed by the linear combination of numerical orthogonal polynomials, regardless of the local surface aperture shape. The entire freeform surface was derived by the overlapping averaging approach. As shown in Fig. 7, a freeform optical element over an irregularly shaped aperture has steep variations at the edge. When approximately fitting it directly by Zernike polynomials over the entire aperture, the local fitting error is relatively large. However, its local deformations at the edge can be characterized finely by the proposed combination method.

## 4.

## Discussion

The benefits of optical freeform surfaces in optical engineering result from the collaborative development of the manufacturing, testing, design tools, and application requirements of freeform surfaces in the past 10 years. Advances in the representation techniques of optical freeform surfaces can further accelerate the improvement of freeform optics in practical applications. From the previous analysis of different types of optical freeform surface representation techniques, we have a large map to help understand the general representation techniques.

Nevertheless, we also need to know how to choose the representation technique of a freeform surface in a specific application. In general cases, when the freeform surface used in the freeform optical design or freeform optical elements is under manufacturing and testing, it can be qualitatively described as a mild or strong freeform surface according to the slope variation degree. For example, Gaussian-like local areas in a freeform surface can be considered a strong freeform surface. The freeform surface shown in Fig. 5 could be considered a mild freeform surface. For a mild freeform surface, diverse representation techniques can be employed according to its practical use. When considering a freeform surface with constrained slopes in its design, Q-type polynomials could be a good choice. For the estimation of a freeform surface with noncircular regularly shaped aperture, orthogonal polynomials with the corresponding aperture shape or numerical orthogonal polynomials are employed. When local surfaces in a few specific fields must be controlled, splines, or radial basis functions may be used. For a freeform surface with large slope variation, hybrid or combination methods could be applied to accurately characterize its local features. On the other hand, in the process of forming the freeform surface (from design to manufacturing to testing), the current technologies of freeform optics also limit the choice of freeform surface representation techniques in practical applications.

Table 5 presents a comparison of several typical freeform surface representation techniques and their advantages. Orthogonal polynomials, because of their orthogonality, completeness, and their aberration-related performance, are frequently applied in freeform surface representation and estimation. We must also note that the orthogonality applies over a continuously defined domain; namely, the sampled data points in the effective aperture must be sufficient. On the other hand, the higher-order polynomial terms can be derived by iterative recurrence. Regarding nonorthogonal polynomials, the XY-type polynomial is always used in freeform surface design, as discussed in Sec. 2, owing to its simplicity compared with orthogonal polynomials. Although the spline function and radial basis function control local features in the freeform surface representation, their ray-tracing efficiency must be improved. These two functions can be integrated into design tools by the “user-defined surface” feature. In the method using a linear combination of base polynomials fitting a strong freeform surface with steep local variation with nanometer accuracy, the use of hundreds or thousands of terms is required. Base polynomials have evident lateral resolving capacity; however, many polynomial terms reduce computational efficiency. In this situation, the hybrid or combination method, discussed in Sec. 3.3, could be an alternative approach to freeform surface representation with finite polynomial terms and high accuracy. Although local fitting coefficients are meaningless for the entire surface, the hybrid or combination method can be applied to characterize local surface features finely for freeform surface estimation.

## Table 5

Comparison of typical freeform surface representation techniques.

Freeform surface | Representation techniques | Advantages | Remarks | |
---|---|---|---|---|

Mild freeform surface | Orthogonal polynomials | Zernike circular polynomials | (1) Each term is related to a specific aberration | (1) The definition domain of the orthogonal polynomial is continuous, not for discrete data points |

Q-type polynomials | (2) Orthogonality and completeness | (2) The derivation of higher-order polynomials is not easy | ||

Zernike polynomials over a regularly noncircular aperture by Gram–Schmidt orthogonalization | (3) One-to-one correspondence with the aperture shape | |||

2-D Chebyshev polynomials | ||||

2-D Legendre polynomials | ||||

Nonorthogonal functions | XY-type polynomials | (1) Related to the aberration | Commonly used in freeform optical design | |

(2) Global definition of freeform surface | ||||

Spline function | (1) Local definition of freeform surface | Ray-tracing with the design tool is challenging | ||

(2) Flexible in controlling surface shapes | ||||

Radial basis function | (1) Local definition of freeform surface | Must be integrated into the design tool using the “user-defined surface” feature | ||

(2) Multicentric feature | ||||

Numerical orthogonal polynomials | (1) Adapt to general aperture shapes | Numerical orthogonal transformation method | ||

(2) Adapt to discrete data points | ||||

Freeform surface with strong slope variation | Linear combination of hundreds or thousands of basis polynomials | Lateral resolving capacity of basis polynomials | Computational efficiency reduced with many polynomial terms | |

Hybrid or combination method | (1) Local feature characterized finely | Local fitting coefficients are meaningless | ||

(2) Only with finite number of polynomial terms |

## 5.

## Conclusion

In the past 10 years, the optical freeform surface representation technique is a hot research topic. The representation techniques of optical freeform surfaces are diverse and not unique in specific applications. The development of representation techniques, the requirements of design and application, and the advancement of manufacturability and testability of optical freeform surfaces will accelerate for the improvement of freeform optics. By focusing on optical freeform surface representation techniques, we review comprehensively the current state of the art of representation techniques and their typical applications. Different types of representation techniques are discussed and compared, which will provide solutions for how to represent an optical freeform surface and which representation technique could be used in a specific application.

Generally, a freeform surface is represented as an explicit mathematical function; it is decomposed as a linear combination of basis function added to the base surface. The base surface is usually a conic surface or sphere. Basis functions can be XY-type polynomials, Zernike polynomials, Q-type polynomials, radial basis functions, or NURBS functions. The differences and advantages of these methods are presented in Table 5. For a strong freeform surface with steep slope variation in local areas, a hybrid or combination method could be applied to accurately characterize local features. In addition, by effectively combining different types of functions, new functions could be constructed in the future for representing optical freeform surfaces. Additionally, the gradient variation of a freeform surface must be considered within the current manufacturing and testing capacities.

Regarding the discrete data points obtained by the direct design method for an optical freeform surface, the obtained discrete data point cloud characterizing the freeform surface must be fitted accurately by an appropriately chosen analytical function to obtain a suitable design starting point for further efficient optimization. A similar process for freeform surface analysis from discrete measured data exists in freeform surface testing. Furthermore, in the design of a freeform optical system with a large field of view, the use of a sectional or stitched freeform surface is a promising technique. The representation technique, as well as the fabrication and measurement, of this type of freeform surface could be challenging.

On the other hand, the functions of most surface types have been integrated into commercial optical design software. However, special functions, such as the radial basis function and the stitched surface, are not included. The “user-defined surface” feature of optical design software is very powerful and flexible and provides more degrees of freedom for freeform optics design. Nevertheless, the problem of the ray-tracing speed and optimization efficiency must be overcome in the future.

## Acknowledgments

The authors would like to thank the Talent Launch Fund of Nanjing University of Information Science and Technology and the National Natural Science Foundation of China (Nos. 61377015, 61505080, and 51705271) for financial support. The authors also want to thank the editor and reviewers for their valuable comments and suggestions.

## References

## Biography

**Jingfei Ye** received his BS degree and PhD in optics from Nanjing University of Science and Technology in 2010 and 2016, respectively. Currently, he is a lecturer in the Department of Optoelectronic Engineering at Nanjing University of Information Science and Technology. His current research interests include optical design and testing especially for freeform optics and optical remote sensing especially for lidar.

**Lu Chen** is currently a PhD candidate at the School of Electronic and Optical Engineering of Nanjing University of Science and Technology. She is working on high-precision metrology and freeform optics.

**Xinhua Li** is a PhD candidate at the School of Electronic and Optical Engineering of Nanjing University of Science and Technology. She is also an associate professor at Jinling Institute of Technology. Her main research is in the field of visual optics.