## Abstract

Many illumination applications require redistributing the irradiance distributions of LED sources with large ray bending. The problem becomes even more challenging for a compact design where the LED size is no longer ignorable. We tackle this problem by simultaneously designing two freeform optical surfaces. An iterative wavefront tailoring (IWT) method is adapted for obtaining the entrance and exit base freeform surfaces with a predefined ray bending regulation under stereographic coordinates (*u*, *v*). The simulated annealing (SA) algorithm is employed for deforming the two base freeform surfaces using the '*uv*' polynomials with the purpose of minimizing the relative root-mean-square deviation (RRMSD) between the simulated irradiance distribution and the prescribed one. The optimizations are implemented in an automated workflow which links the optimization engine, 3D modeling software and ray tracing software. The effectiveness of the proposed method is illustrated by designing several double-freeform-surface lenses (central heights: 10 mm) with different ray bending regulated base surfaces and 10-*th* order *uv* polynomial departures for generating 500 × 200 mm^{2} uniform irradiance distributions at a distance of 100 mm from 2 × 2 mm^{2} and 3 × 3 mm^{2} sources, respectively.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Typical LED sources have planar Lambertian emitting surfaces that emit light hemispherically and generate very non-uniform illumination patterns on the targets perpendicular to the optical axes. Secondary reflective or refractive optics are usually employed to transform the LED light distributions into prescribed ones which satisfy certain lighting application requirements. Here we focus on secondary lenses which are more commonly used. For a rotationally symmetric design, a 'mushroom' lens can be employed to uniformly spread the LED light distribution [1]. Freeform lens is necessary when the irradiance distribution shape of the target is different from that of the LED source. A large ray bending ability is required for the freeform lens when the deviation between the source and target irradiance distributions is very large (see Fig. 1). Difficulty also arises from the fact that the increasing applications of LEDs in more and more compact systems make the source sizes non-negligible.

Compared with the most commonly used freeform lenses with single freeform surfaces, freeform lenses with double freeform surfaces, i.e., both the entrance and exit surfaces are freeform, are more powerful for bending the light rays with large deflection angles. Parkyn [2,3] probably pioneered the use of double freeform surfaces in the cases of large ray bending, where the intermediate ray vectors between the entrance and exit freeform surfaces are set to be proportional to the average of the input and output ray vectors. Bäuerle *et al*. [4] proposed a least squares construction of double freeform surfaces following an approximate ray mapping computed by the optimal transport. Bruneton *et al*. [5] proposed a ray mapping optimization method in the double freeform surfaces design considering the light flux difference, surface smoothness and distance to target’s boundary as the merit function. Wu *et al*. [6] did a direct formulation of double freeform surfaces, where the entrance freeform surface is predefined for producing a pattern of which the distribution range is closed to that of the target irradiance distribution. Instead of a direct determination, Feng *et al*. [7] proposed an iterative wavefront tailoring (IWT) method to simplify the freeform lens design through the intermediate construction of an outgoing wavefront. This method is capable of constructing a single freeform surface lens with two freeform surfaces where the entrance freeform surface is predefined by an analytical expression. All these design methods of double freeform surfaces are developed under point source approximations, which will become less valid when the extended sizes of the light sources must be taken into account.

There have been several methods devised for extended light sources. Over compensation methods iteratively apply over compensations to the point source designs [8–12]. Deconvolution methods reduce the blur effects caused by extended light sources applied on the point source designs [13–16]. The Simultaneous Multiple Surface (SMS) method can design three freeform surfaces to couple three edge wavefronts from the three corners of a square LED chip to the corresponding three wavefronts of the target intensity pattern [17]. This wavefront coupling scheme is integrated into a wavefront tailoring technique, which can design a highly compact double-freeform-surface lens delivering a flat-top illumination pattern [18]. The edge ray mapping method can design a single freeform surface lens for generating a convex polygonal illumination pattern from an extended light source based on a mapping of the maximally off-axis rays from the the edges of the source through the edge of the lens to the edges of the target polygon [19]. Energy accumulating optimization based on dynamic programming is proposed to design a single discrete freeform surface that is composed of sufficiently small pierces [20]. Another interesting work is the employment of the expectation maximization algorithm for arranging the pinhole images of the extended source generated by points on the optical surface to form a prescribed irradiance pattern [21].

A very simple and versatile way to solve the irradiance tailoring problem with an extended source is to represent the freeform surfaces using basis functions and optimize their weights to minimize the deviation between the simulated irradiance distribution and the prescribed one (see, e.g., [22–29]). In fact, this parametric optimization way has dominated imaging optics design for a long time. For illumination optics design, parametric optimization is gaining more and more attention as the improvements in computing power and developments of highly efficient ray tracing techniques. One of the key issues in parametric optimization is to find an appropriate surface description technique. The simplest analytic expression for a freeform surface might be the $XY$ polynomials [22], which, for instance, have been applied in designing a single freeform surface [23] and double freeform surfaces [24] respectively for producing rectangular illumination patterns from extended LED sources. The drawback of expressing a freeform surface with a $XY$ polynomial is that the light rays emitted from a typical LED source into semi-space cannot be fully captured. Freeform surface representation techniques that can cover the hemispherical light emissions of LED sources include Spherical harmonics [25], B-splines [26–28] and T-spines [29]. Another interesting surface representation is based on the superposition of several base surfaces generated from pedal curves, which is integrated with a relative illumination calculation using étendue to evaluate tolerance [30,31]. Very few of theses studies involve tackling large ray bending cases. Recently, Byzov *et al*. [28] demonstrate an impressive example that double B-spline surfaces can largely bend the light rays of an extended LED source to generate a square uniform illumination pattern, where a local optimization algorithm implemented in the *fminunc* MATLAB routine is employed to optimize the B-spline parameters and a self-developed ray tracing technique is used to calculate the merit function.

Here, we propose to optimize double freeform surfaces simultaneously using simple polynomials with respect to stereographic coordinates $(u,v)$. Instead of directly describing the complete double freeform surfaces, the $uv$ polynomials are employed to describe the radius distance deviations from two base freeform surfaces. The IWT method is adapted for designing the entrance and exit base freeform surfaces with a predefined ray bending regulation. The SA algorithm, one of the heuristic optimization methods, is chosen as the optimizer to find a global minimum of merit function with a larger probability than the local optimization methods. An automated workflow that combines optimization engine, 3D modeling software, and ray tracing software is provided to implement the optimization. This method is detailedly described in Section 2. In Section 3, we provide freeform lens designs which are aimed for generating $500\times 200$ mm$^2$ uniform irradiance distributions at a distance of 100 mm from $2\times 2$ mm$^2$ and $3\times 3$ mm$^2$ sources respectively. The influences of the base freeform surfaces with different ray bending regulation parameters, the number of polynomial terms and the ray tracing settings on the optimizations are also discussed. Finally, the conclusion and a short discussion are provided in Section 4.

## 2. Design method

#### 2.1 General statement

Our goal is to design a single lens with double freeform surfaces that can produce a desired irradiance distribution on a given target plane at $z=H$ from an extended light source (see Fig. 2(a)). The prescribed irradiance distribution is denoted by $E_t(\xi , \eta )$, where $(\xi , \eta ) \in \Gamma \subset \mathbb {R}^2$. The light source is assumed to have some arbitrary intensity distributed in the semi-space, and thus the Cartesian coordinates are not suitable here. $XY$, Zernike and other polynomials for Cartesian coordinates are no longer applicable for describing freeform optical surfaces that cover the hemispherical light emissions. Instead of spherical coordinates, we prefer using the stereographic projection which establishes a correspondence between a point $(X, Y, Z)$ on the unit sphere and a point $(u,v)$ on the plane $Z=0$ based on the projection from the south pole $(0, 0, -1)$ (see Fig. 2(b)). $(X, Y, Z)$ can be described using $(u,v)$ as:

The stereographic coordinates $(u,v)$ belong to a unit circular domain $\Sigma =\{(u,v)|u^2+v^2\leq 1 \}$ which corresponds to the upper hemisphere. As shown in Fig. 2(c), the $(u,v)$ coordinates in the unit circular domain $\Sigma$ can be transformed from the coordinates $(u_0,v_0)$ in a square domain $\tilde {\Sigma }=\{(\tilde {u},\tilde {v})|-1 \leq \tilde {u},\tilde {v}\leq 1\}$ through:

We want to keep the form of $XY$ polynomials which is simple and efficient for optimization. However, unlike a $XY$ polynomial which describes the height function $z$, we introduce a composite representation to describe the radius distance function of a freeform surface that can cover the hemispherical light emission:

For our design problem with double freeform surfaces, Eq. (3) can be used to represent the radial distance function of the entrance freeform surface. Similarly, the radial distance function of the exit freeform surface can be described as:

where $\rho _0'(u,v)$ refers to the radius distance function of the exit base surface. The position vectors $(x',y',z')$ on the exit surface can also be obtained from Eq. (4) by replacing $\rho (u,v)$ with $\rho '(u,v)$.For a case of large ray bending, the required freeform surfaces are probably be severely deformed. It could be helpful if we can provide some base freeform surfaces that are close to the required freeform surfaces. In the following, we will first introduce the construction of the base freeform surfaces using a point source design method and then describe the optimization strategies associated with an automated workflow.

#### 2.2 Base surfaces design

Here we introduce the design of the entrance and exit base surfaces based on an IWT method under stereographic coordinates [32] integrated with a double-surface reconstruction procedure [12].

Suppose the approximated point source has an intensity $I(u,v)$, then the stereographically projected irradiance (SPI) $E(u,v)$ on the $(u,v)$ plane can be obtained as [33]:

We implement the IWT procedure under the Cartesian coordinates $(\tilde {u},\tilde {v})$. The SPI for $(\tilde {u},\tilde {v})$ can be obtained by coordinates transformation [32]:

We first apply energy conservation between source and target in differential form as:

After that, we will make a connection between Eq. (9) and the gradients of an outgoing wavefront $\mathbf {W}=(s, t, w)$ that is next to the exit surface. According to Fermat’s principle and chain’s rule, the gradients of $w$ with respect to $\tilde {u}$ and $\tilde {v}$ can be expressed as [32]:

Form this equation, we express the target coordinates $(\xi , \eta )$ as [32]:

We then insert Eq. (11) into Eq. (9) to obtain a MA equation of $w(\tilde {u} ,\tilde {v})$ [32]

Next, we will describe the reconstruction of double freeform surfaces denoted by position vectors $\mathbf {P}_{i,j}$ and $\mathbf {Q}_{i,j}$ following the ray map based on [12]. After discretizing $\tilde {\Sigma }$, the grid points $(u_{i,j},v_{i,j})$ on the discretized domain $\Sigma$ can be immediately generated from Eq. (2). Then, an incident ray sequence can be obtained as:

After specifying an incident ray sequence, the double-surface reconstruction can be more convenient if we could know the corresponding outgoing ray sequence. However, the outgoing ray sequence is coupled with the surface data. Therefore, we employ an iterative procedure to simplify the surface reconstruction.

We first give an estimate of the outgoing ray sequence as:

**Unit**refers to the operation of obtaining the unit vector. This estimate is based on the far field approximation that the lens size can be omitted. We then define a unit intermediate ray sequence $\hat {\mathbf {R}}_{i,j}$ within the lens as [3]:

The physical meaning of Eq. (17) is that a chord that links two adjacent surface points is perpendicular to the average of the normals at the two points [33]. The surface points $\mathbf {P}_{i,j}$ can be represented as $\rho _{0,i,j}\hat {\mathbf {I}}_{i,j}$ which are inserted into Eq. (17) to obtain a system of linear equations of $\rho _{0,i,j}$ [33]:

Equation (18) can be rewritten as a single matrix equation, and Hermann’s method [35] can be used to obtain a least squares solution. After obtaining an entrance surface based on least squares, the interior ray sequence $\hat {\textbf {R}}_{i,j}$ is recalculated by first computing the normal vectors to the entrance surface and then applying Snell’s Law. We then compute the required normal vectors $\hat {\textbf {N}}_{i,j}'$ of the exit surface based on Snell’s law:

The exit surface points $\mathbf {Q}_{i,j}$ can also be obtained based on least squares, where the basic relationships are:

Similarly, Eq. (20) can be rewritten as a linear system of equations of the distances between $\mathbf {P}_{i,j}$ and $\mathbf {Q}_{i,j}$, which can also be solved using Hermann’s method [35]. After specifying an exit surface $\mathbf {Q}_{i,j}$, the output ray vectors can be updated as:

The updated output ray vectors are inserted into Eq. (15) to repeat the double-surface reconstruction until a stop criterion is reached.After the reconstruction of the entrance and exit surfaces following the above procedure, we can reconstruct an intermediate wavefront $\mathbf {W}_{i,j}=(s_{i,j}, t_{i,j}, \tilde {w}_{i,j})$ using a least squares method, where $\tilde {w}$ is used to differentiate from $w$ in Eq. (12). The basic relationships to implement the least squares are:

Equation (22) can be rewritten as a linear system of equations of the distances between $\mathbf {Q}_{i,j}$ and $\mathbf {W}_{i,j}$, and Hermann’s method [35] is employed again to acquire a least squares solution. The updated $(s_{i,j}, t_{i,j})$ are inserted into Eq. (12) to solve a new wavefront $w_{i,j}$ and thus obtain a new ray map $(\xi _{i,j}, \eta _{i,j})$ based on Eq. (11). The new ray map can help reconstruct a new entrance surface and exit surface. We iteratively revise the ray map based on solving Eq. (12) and accordingly update the double freeform surfaces using the iterative reconstruction procedure to improve performance. A multi-scale strategy could be used to speed up the computation [7].

It is necessary to note that the source domain for the IWT computation is usually set as $\tilde {\Sigma }_\gamma =\{(\tilde {u},\tilde {v})|-\gamma \leq \tilde {u},\tilde {v} \leq \gamma \}$, where $\gamma$ is a positive value smaller than 1. Since $\tilde {\Sigma }_\gamma$ is smaller than $\tilde {\Sigma }=\{(\tilde {u},\tilde {v})|-1 \leq \tilde {u},\tilde {v} \leq 1 \}$, its corresponding domain $\Sigma _\gamma$ after coordinates transformation based on Eq. (2) becomes smaller than the unit circular domain $\Sigma$ of which the boundary corresponds to the incident light rays which are parallel to the $z=0$ plane. Such a strategy could omit the light rays that are almost parallel to the $z=0$ plane, which may experience total internal reflection and result in a failure of the IWT computation [32]. A rigorous determination of $\gamma$ could be very difficult. We take a trial and error strategy which gradually scales down the $\gamma$ value to check if a reasonable solution could be obtained.

The computed entrance and exit surfaces can no longer cover the hemispherical emission from a typical LED source. However, an initial entrance and exit surfaces that can cover the hemispherical emission could be beneficial for the subsequent surface representation and optimization. For the entrance surface, the computed $\rho _0$ values for the discretized domain $\Sigma _\gamma$ are interpolated and extrapolated to those values for the discretized unit circular domain $\Sigma$ using the Matlab *scatteredInterpolant* function. For the computed exit surface, we first obtain the $\rho '_0$ values as: $\rho '_{0,i,j}=(x_{0,i,j}'^2+y_{0,i,j}'^2+z_{0,i,j}'^2)^{1/2}$, where $x_{0,i,j}$, $y_{0,i,j}$ and $z_{0,i,j}$ are the three coordinates of $\mathbf {Q}_{i,j}$. The corresponding stereographic coordinates $(u'_{i,j}, v'_{i,j})$ can be obtained as:

The $\rho '_0$ values with respect to $(u'_{i,j}, v'_{i,j})$ are interpolated and extrapolated to those values for $(u_{i,j}, v_{i,j})$, where $(u_{i,j}, v_{i,j})$ are the coordinates of the grid points of the discretized domain $\Sigma$. Finally, we obtain an entrance base surface and exit base surface both of which can cover the hemispherical emission from an LED source.

#### 2.3 Optimization

After obtaining the entrance and exit base surfaces, we discuss the optimization strategies of the $uv$ polynomials. The optimization workflow is shown in Fig. 3. The optimizer can generate certain values of the weights $\mathbf {a}=[a_{i,j}]$ and $\mathbf {b}=[b_{i,j}]$, which are inserted into Eq. (3) and Eq. (5) to obtain the radius lengths $\rho _{i,j}$ and $\rho '_{i,j}$, respectively. The scattered points on the entrance and exit freeform surfaces, respectively, can be specified based on Eq. (4). The scattered data is inserted into a 3D moldeling software to obtain a solid model of the freeform lens. This freeform lens model is imported into a ray tracing simulation software. The simulated irradiance distribution is used to calculate the merit function (MF) value. With the MF value, the optimizer updates the values of the weights and restarts a new iteration. The optimizer also stops the optimization after some stop criteria are satisfied. In the following, we provide some details for implementing the optimization.

Since Monte-Carlo ray tracing for irradiance analysis can produce uncertainty or error with gradient computation [36], optimizers employing derivatives are not preferred here. Simplex optimization method [36] could be appropriate for our case, but we choose the SA method because it may find the global minimum with a larger probability. SA is a random searching optimization algorithm which is developed based on the similarity between the annealing process of solid matter and the general combinatorial optimization problem [37]. The SA algorithm starts from an initial temperature $T_0$. With the continuous decrease of temperature $T$ and the characteristic of the probability jumping, the SA algorithm searches for the global optimal solution of the MF randomly in the solution space. Here, the SA algorithm is performed in Matlab using the *simulannealbnd* function. The probability of accepting a worse result is $1 \big /\left [1+\textrm {exp}\left (\Delta / \textrm {max}(T)\right )\right ]$, where $\Delta$ denotes the difference between the new and old MF values [38]. A well selected initial temperature $T_0$ is beneficial. If $T_0$ is too large, meaningless jumps of variables when the optimizer searching the solution space could consume a large amount of time. If $T_0$ is too small, the optimizer may be trapped into local optimality easily.

We employ the relative root-mean-square deviation (RRMSD) between the simulated and target irradiance distributions as the MF:

where $\mathbf {E}_t$ and $\mathbf {E}_s$ denote the matrices of the target and simulated irradiance distributions respectively, and $\Vert \cdot \Vert _F$ is the Frobenius norm. Besides, direct controls on geometric relations, e.g., the radius distances of the second surface should be larger than those of the first surface and the bottom of inner surface doesn’t intersect with the source, can be implemented through applying penalties to the MF.To realize an automated optimization, we link three types of software. Matlab is used as the execution control platform and optimization engine. Rhinoceros is employed for obtaining a solid lens model from the scattered points data generated from Matlab. LightTools is used for irradiance analysis based on Monte-Carlo ray tracing. Matlab can read the data of the simulated irradiance distribution and calculate the MF value. Matlab controls Rhinoceros and LightTools via application programming interfaces (APIs). We need to mention that the surface modeling in Rhinoceros is based on NURBS (Non-Uniform Rational Basis Splines). However, the number of the scattered surface points outputted from Matlab is sufficiently large so that the surface representation inconsistency has a small impact on the simulation results.

## 3. Design examples

In this section, we demonstrate the feasibility of the proposed design method by generating $500 \times 200$ mm$^2$ uniform illumination patterns at $z=100$ mm from $2 \times 2$ mm$^2$ and $3 \times 3$ mm$^2$ Lambertian sources respectively. The ratio of the length, width and height of the illumination reaches 5:2:1, which requires large ray bending abilities for the freeform lenses. The central heights of the entrance and exit base surfaces are set as 3.5 mm and 10 mm respectively. During optimization, the central heights of exit freeform surfaces are fixed at 10 mm, while there are no restrictions on the central heights of the entrance freeform surfaces. The material of the lenses are set as PMMA (Polymethyl methacrylate). We choose the 10th order $uv$-polynomials for deforming both the entrance and exit surfaces. Since the prescribed illumination pattern is axisymmetric, we only need to optimize the even order terms (and the constant term) of the $uv$-polynomials. Therefore, there are 42 variables to be optimized (21 for each), and their initial values are all set to zero. After connecting Matlab, Rhinoceros and LightTools, the SA optimizer automatically deforms the entrance and exit surfaces by altering the 42 variables to check the variation of the MF value until satisfying a certain number of iterations which is set as 2,000 here. The number of the $(u,v)$ grid points is set as $101 \times 101$. Therefore, in each MF evaluation, the number of the outputted scattered points of either the entrance or the exit surface is $101 \times 101$, which is sufficiently large so that the surface modeling using NURBS (in Rhinoceros) has a very small deviation from that using $uv$-polynomials. For the simulation step, 2 million rays are traced for each MF evaluation, and a resolution of $100 \times 100$ for the simulated irradiance distributions within the target domain is set to compute the MF value. Here the MF is set to be 10 times of the RRMSD.

The designs for the $2 \times 2$ mm$^2$ and $3 \times 3$ mm$^2$ sources share the same base surfaces which are calculated by the IWT method in a multi-scale way. The IWT algorithm is first implemented in a very coarse grid of $16 \times 16$. A Cartesian grid of the target domain $\Gamma =\{(\xi , \eta )| -100\ \textrm {mm} \leq \xi \leq 100\ \textrm {mm}, -250\ \textrm {mm} \leq \eta \leq 250\ \textrm {mm}\}$ with equal spacings along $x$ and $y$ axes respectively is set as the initial map. This initial map can be used to acquire an estimate of the outgoing wavefront that is next to the exit base surface. The outgoing wavefront estimate can be inserted into the wavefront equation shown in Eq. (12) to obtain a updated ray map. The computed ray map after two iterations is then interpolated into the size of $32 \times 32$, which is used to start the algorithm on the $32 \times 32$ grid. Such a process is repeated until finishing the computation on the $256 \times 256$ grid. In accordance with the discretization of the $uv$ polynomial surfaces, the computed $\rho _0$ and $\rho _0'$ values are interpolated and extrapolated into the $(u,v)$ grid of $101 \times 101$.

Figures 4(a)–4(d) show the initial double-freeform-surface lenses designed by the mutlti-scale IWT algorithm with different ray bending regulations of $\alpha$=0.4, 0.5, 0.6 and 0.7. Figure 4(e) also provides a single freeform surface lens ($\alpha$=1), where the inner surface is a hemisphere with the radius of 3.5 mm. We can see that the lens becomes more and more compact as $\alpha$ decreases i.e. the inner surface plays a more and more important role in light deflection. The source domains involved in the IWT computations for $\alpha$=0.4, 0.5 and 0.6 are all set as $\tilde {\Sigma }_{\gamma =0.9}=\{(\tilde {u},\tilde {v})|-0.9 \leq \tilde {u},\tilde {v} \leq 0.9 \}$. When $\alpha$ is set as 0.7, we employ a slightly smaller source domain $\tilde {\Sigma }_{\gamma =0.86}$ for omitting more light rays with large polar angles which may experience total internal reflection. Although the source domain for $\alpha$=0.4, 0.5, 0.6 or 0.7 has a reduction of more than $10\%$ from the domain $\tilde {\Sigma }$, only very small amount of flux is not included in the calculation since the light source has a Lambertian intensity. For the single freeform surface design of $\alpha =1$, we shrink the source domain much further to $\tilde {\Sigma }_{\gamma =0.68}$ for acquiring reasonable results.

Figures 4(f)-(j) provide the simulated irradiance distributions of the five initial freeform lenses shown in Figs. 4(a)–4(e) respectively, with a point-like source (size: $10^{-2} \times 10^{-2}$ mm$^2$). All these irradiance distributions look uniform and sharp at the edges, and their RRMSD values are less than $0.14$ (see Table 1). We can also see from Table 1 that each double-freeform-surface lens performs a high light efficiency (the ratio of the flux projected into the target domain $\Gamma$ to the total flux of the source). The light efficiency of the single freeform surface lens for $\alpha =1$ is almost ten percentage points lower than each double-freeform-surface lens when considering Fresnel losses.

Figures 4(k)–4(o) show the simulated irradiance distributions of the five point-source designs when we replace the point-like source with a $2 \times 2$ mm$^2$ one. We can clearly observe the deformations influenced by the extended size of the source on the point-source designs: hot regions show up (especially for $\alpha$ = 0.4 and 0.5) and the edges become less sharp. We can also observe from Figs. 4(k)-(o) that the extended size of the light source has a larger influence on the irradiance profiles along $y$ directions than those along $x$ directions. The RRMSD values are shown in Table 2 (Starting RRMSD), and we can see that the point source design for $\alpha =1$ has the largest RRMSD value of 0.2301. However, the point source design for $\alpha =0.4$ has the biggest increase in the RRMSD value, from 0.1146 to 0.2256. The performances are dropped further for the five point-source designs when we increase the source size to $3 \times 3$ mm$^2$, as shown in Figs. 4(p)–4(t) and Table 2. The point source design for $\alpha =0.4$ is most affected and its RRMSD value is as high as 0.3415. The point source design for $\alpha =1$ experiences a smallest increase in the RRMSD value when the $2 \times 2$ mm$^2$ source is changed into a $3 \times 3$ mm$^2$ one, which is probably a result that the freeform lens for $\alpha =1$ (see Fig. 4(e)) has the biggest size.

Figures 5(a)–5(e) show the optimized freeform lenses for the $2 \times 2$ mm$^2$ source with different base surfaces (see Figs. 4(a)–4(e)) of $\alpha$=0.4, 0.5, 0.6, 0.7 and 1. The inner surface for the case of $\alpha$=1 is kept as a hemisphere during optimization. Figures 5(f)–5(j) show the corresponding simulated irradiance distributions, where we can clearly observe the performance improvements compared with Figs. 4(k)–4(o). The changes of the RRMSD values and the light efficiencies before and after optimizations are shown in Table 2. We can see that the RRMSD values for the double-freeform-surface designs have been effectively reduced to be lower than 0.1307. The light efficiency for each design has a little reduction because a small amount of the flux leaks out of the target domain. Compared with double-freeform-surface designs, the single freeform surface design has the worst RRMSD value. In addition, the light efficiency of the single freeform surface design is only around $75\%$ which is about 9 percentage points lower than that of each double-freeform-surface design. These two facts indicate that a conventional design with an entrance hemisphere and an exit freeform surface is probably be insufficient for the large ray bending case here. The optimized freeform lenses for the $3 \times 3$ mm$^2$ source are shown in Figs. 5(k)–5(o), and their corresponding simulated irradiance distributions are shown in Figs. 5(p)–5(t). The RRMSD values for the double-freeform-surface designs still experienced effective reductions (see Table 2). Although the initial double freeform lens for the case of $\alpha$=0.4 has the largest RRMSD value, its corresponding optimized double freeform lens has a better performance than the optimized single freeform lens. Compared with the results of the $2 \times 2$ mm$^2$ source, the reductions of the RRMSD values and light efficiencies for the the $3 \times 3$ mm$^2$ source are consistent with the intuition that it is more difficult to handle the irradiance tailoring problem for a larger extended source.

Figure 6 shows the $y=0$ and $x=0$ cross-sectional profiles of the base freeform surfaces and the optimized freeform surfaces with the $2 \times 2$ mm$^2$ and $3 \times 3$ mm$^2$ sources for the five cases of $\alpha$=0.4, 0.5, 0.6, 0.7 and 1. For each case, it seems that the inner surface is higher for a larger source when the central height of the exit surface maintains at 10 mm. It can also be seen from Fig. 6 that the $x=0$ profiles have bigger changes than the $y=0$ profiles as $\alpha$ varies. The $x=0$ profiles of the exit surfaces are convex for $\alpha$=0.4, but their inner sections become more and more concave as $\alpha$ increases. The barb-like structure at the bottom of the inner surface profile shown in Fig. 6(d) or 6(f) adds difficulties to the lens fabrication. However, this structure can be modified to facilitate fabrication, and such a modification has a small influence on the optical performance because the corresponding part receives a very small amount of the source flux. A fundamental solution to this problem is to add the manufacturability as a constraint in the optimization.

Take the case of $\alpha =0.5$ and $3 \times 3$ mm$^2$ source for example, we analyse the influence of different orders of the $uv$ polynomials on the optimization performances. Figure 7 show the $y=0$ and $x=0$ cross-sectional profiles of the optimized freeform surfaces that include 8, 10 and 12 orders $uv$ polynomials respectively, and their differences are almost invisible. The resulting RRMSD values and light efficiencies are very close (see Table 3). The reason is probably that the point source designs provide good base surfaces.

We also provide in Table 4 that the influence of the number of rays and receiver bins on the optimization behaviour ($\alpha =0.5$; source size: $3 \times 3$ mm$^2$). As can be seen from Table 4 that, the RRMSD values and light efficiencies experience small changes when the number of rays is altered into $1 \times 10^6$ or $4 \times 10^6$ while the number of bins maintains at $100 \times 100$. The error estimate at the peak is decreased from $11.25\%$ to $6.14\%$ when we increase the number of rays from $1 \times 10^6$ to $4 \times 10^6$. This is in accordance with the Rose model that the signal-to-noise ratio (SNR) is proportional to the root of the number of rays [39]. Table 4 also illustrates that the change of the receiver bins (grid size of the target domain) has a little effect on the RRMSD and light efficiency but has a big effect on the error estimate. As observed in Table 4, when we reduce the number of bins while maintaining the number of rays, the error estimate at the peak can also be reduced. This is because that each bin will collect more rays when the number of the bins is decreased. We need to mention that mesh smoothing with an integration width of 3 bins is applied for each irradiance evaluation in LightTools. Such an operation can reduce the sampling noise that is inherent in the Monte-Carlo ray tracing.

## 4. Conclusion

We have proposed a double-freeform-surface lens design method for tackling the large ray bending irradiance tailoring problem with an extended (LED) source. The entrance and exit freeform surfaces are specified with base freeform surfaces for point source and additional $uv$ polynomial departures. The base freeform surfaces are designed by an IWT method with a predefined ray bending regulation between the entrance and exit surfaces. The SA algorithm is employed for optimizing the $uv$ polynomial deformations in an automated workflow that joins the optimization engine, 3D modeling software and ray tracing software. With different regulated base freeform surfaces and 10 order $uv$ polynomial departures, we provide challenging designs aiming at generating $500 \times 200$ mm$^2$ uniform irradiance distributions at $z$ = 100 mm. Results show that we can acquire double-freeform-surface lenses which produce illumination patterns with RRMSD values less than $0.131$ and light efficiencies higher than $91.9\%$ (without Fresnel losses) for a $2 \times 2$ mm$^2$ Lambertian sources, and acquire double-freeform-surface lenses which produce illumination patterns with RRMSD values less than $0.197$ and light efficiencies higher than $87.6\%$ for a $3 \times 3$ mm$^2$ Lambertian source.

Our design method can be extended to design freeform lenses which can generate non-uniform irradiance distributions e.g., super-Gaussian and cosine distributions, and may also be applicable for non-rectangular and non-axisymmetrical irradiance distributions. For non-axisymmetrical cases, odd terms of the $uv$-polynomials should be involved in the optimization. One of the shortcomings of our method is that the optimization is slow. The total computation time for each optimization could take more than 16 hours on an Intel Core i9-10900k at 3.7 GHz with 128 GB RAM. Ray tracing spends most of the time. When the number of rays was changed from $2 \times 10^6$ to $4 \times 10^6$ (see Table 4), more than 25 hours were costed to finish the optimization. Connections between Matlab, Rhinoceros and LightTools also take up a certain amount of time. Accelerated ray tracing mode in LightTools can be selected to balance speed and accuracy. Hirst *et al.* developed a deterministic method for irradiance evaluation from an extended Lambertian source through a smooth single refractive optical surface [40]. This method could be extended to tackle the irradiance computation problem of a double-freeform-surface lens. Future work may also include the incorporation of manufacturability and tolerance issues into the optimization process, the applications of orthogonal polynomials with respect to $(u,v)$ and the generation of complicated illumination patterns.

## Funding

National Key Research and Development Program of China (2017YFA0701200); National Natural Science Foundation of China (11704030).

## Acknowledgment

The authors are grateful for the educational licenses of Matlab, Rhinoceros and LightTools. Zexin Feng thanks Xianglong Mao for his help in linking the three types of software.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **W. A. Parkyn and D. G. Pelka, “Auxiliary lens to modify the output flux distribution of a TIR lens,” U.S. patent 5577493 (1996).

**2. **W. A. Parkyn, “Illumination lenses designed by extrinsic differential geometry,” Proc. SPIE **3482**, 389–396 (2011). [CrossRef]

**3. **B. Parkyn and D. Pelka, “Free-form lenses designed by a pseudo-rectangular lawnmower algorithm,” Proc. SPIE **6338**, 633808 (2006). [CrossRef]

**4. **A. Bäuerle, A. Bruneton, R. Wester, J. Stollenwerk, and P. Loosen, “Algorithm for irradiance tailoring using multiple freeform optical surfaces,” Opt. Express **20**(13), 14477–14485 (2012). [CrossRef]

**5. **A. Bruneton, A. Bäuerle, R. Wester, J. Stollenwerk, and P. Loosen, “High resolution irradiance tailoring using multiple freeform surfaces,” Opt. Express **21**(9), 10563–10571 (2013). [CrossRef]

**6. **R. Wu, S. Chang, Z. Zheng, L. Zhao, and X. Liu, “Formulating the design of two freeform lens surfaces for point-like light sources,” Opt. Lett. **43**(7), 1619–1622 (2018). [CrossRef]

**7. **Z. Feng, D. Cheng, and Y. Wang, “Iterative wavefront tailoring to simplify freeform optical design for prescribed irradiance,” Opt. Lett. **44**(9), 2274–2277 (2019). [CrossRef]

**8. **Y. Luo, Z. Feng, Y. Han, and H. Li, “Design of compact and smooth free-form optical system with uniform illuminance for LED source,” Opt. Express **18**(9), 9055–9063 (2010). [CrossRef]

**9. **Z. Feng. Y. Luo and Y. Han, “Design of LED freeform optical system for road lighting with high luminance/illuminance ratio,” Opt. Express **18**(21), 22020–22031 (2010). [CrossRef]

**10. **X. Mao, H. Li, Y. Han, and Y. Luo, “Two-step design method for highly compact three-dimensional freeform optical system for LED surface light source,” Opt. Express **22**(S6), A1491–A1506 (2014). [CrossRef]

**11. **R. Wester, G. Müller, A. Völl, M. Berens, J. Stollenwerk, and P. Loosen, “Designing optical free-form surfaces for extended sources,” Opt. Express **22**(S2), A552–A560 (2014). [CrossRef]

**12. **Z. Feng, D. Cheng, and Y. Wang, “Freeform road lighting lens design,” Proc. SPIE **10693**, 106930E (2018). [CrossRef]

**13. **T. Kiser and M. Pauly, “Caustic Art,” EPFL Technical Report (2012).

**14. **D. Ma, Z. Feng, and R. Liang, “Deconvolution method in designing freeform lens array for structured light illumination,” Appl. Opt. **54**(5), 1114–1117 (2015). [CrossRef]

**15. **M. Brand and D. A. Birch, “Freeform irradiance tailoring for light fields,” Opt. Express **27**(12), A611–A619 (2019). [CrossRef]

**16. **C. Bösel and H. Gross, “Compact freeform illumination system design for pattern generation with extended light sources,” Appl. Opt. **58**(10), 2713–2724 (2019). [CrossRef]

**17. **S. Sorgato, R. Mohedano, J. Chaves, M. Hernández, J. Blen, D. Grabovičkić, P. Benítez, J. Miñano, H. Thienpont, and F. Duerr, “Compact illumination optic with three freeform surfaces for improved beam control,” Opt. Express **25**(24), 29627–29641 (2017). [CrossRef]

**18. **S. Sorgato, J. Chaves, H. Thienpont, and F. Duerr, “Design of illumination optics with extended sources based on wavefront tailoring,” Optica **6**(8), 966–971 (2019). [CrossRef]

**19. **D. A. Birch and M. Brand, “Design of freeforms to uniformly illuminate polygonal targets from extended sources via edge ray mapping,” Appl. Opt. **59**(22), 6490–6496 (2020). [CrossRef]

**20. **K. Wang, Y. Han, H. Li, Y. Luo, C. Sun, Z. Hao, B. Xiong, J. Wang, and L. Wang, “Design of high-compactness freeform optical surfaces via energy accumulating optimization,” Opt. Express **24**(26), A1489–A1504 (2016). [CrossRef]

**21. **A. Völl, M. Berens, R. Wester, P. Buske, J. Stollenwerk, and P. Loosen, “Freeform optics design for extended sources in paraxial approximation exploiting the expectation maximization algorithm,” Opt. Express **28**(24), 37004–37014 (2020). [CrossRef]

**22. **P. Benítez and J. C. Miñano, “The future of illumination design,” Opt. Photonics News **18**(5), 20–25 (2007). [CrossRef]

**23. **Z. Liu, P. Liu, and F. Yu, “Parametric optimization method for the design of high-efficiency free-form illumination system with a LED source,” Chin. Opt. Lett. **10**(11), 112201 (2012). [CrossRef]

**24. **J. Lai, X. Li, and P. Ge, “Designing double freeform surfaces for point source and extended sources using a construction, iteration and optimization process,” Lighting Res. Technol. **51**(7), 1118–1127 (2019). [CrossRef]

**25. **C. Gannon and R. Liang, “Using spherical harmonics to describe large-angle freeform lenses,” Appl. Opt. **57**(28), 8143–8147 (2018). [CrossRef]

**26. **Z. Zhuang, P. Surman, and S. Thibault, “Multi-element direct design using a freeform surface for a compact illumination system,” Appl. Opt. **56**(32), 9090–9097 (2017). [CrossRef]

**27. **M. A. Moiseev, S. V. Kravchenko, and L. L. Doskolovich, “Design of efficient LED optics with two free-form surfaces,” Opt. Express **22**(S7), A1926–A1935 (2014). [CrossRef]

**28. **E. V. Byzov, S. V. Kravchenko, M. A. Moiseev, E. A. Bezus, and L. L. Doskolovich, “Optimization method for designing double-surface refractive optical elements for an extended light source,” Opt. Express **28**(17), 24431–24443 (2020). [CrossRef]

**29. **A. S. Isaac and C. Neumann, “Optimization of freeform surfaces using intelligent deformation techniques for LED applications,” Adv. Opt. Technol. **7**(1-2), 67–80 (2018). [CrossRef]

**30. **J. Sasián, D. Reshidko, and C.-L. Li, “Aspheric/freeform optical surface description for controlling illumination from point-like light sources,” Opt. Eng. **55**(11), 115104 (2016). [CrossRef]

**31. **J. Ryu and J. Sasian, “Tolerancing a lens for LED uniform illumination,” Proc. SPIE **10377**, 1037703 (2017). [CrossRef]

**32. **Z. Feng, D. Cheng, and Y. Wang, “Iterative freeform lens design for prescribed irradiance on curved target,” Opto-Electron. Adv. **3**(7), 200010 (2020). [CrossRef]

**33. **Z. Feng, B. D. Froese, and R. Liang, “Freeform illumination optics construction following an optimal transport map,” Appl. Opt. **55**(16), 4301–4306 (2016). [CrossRef]

**34. **R. A. Hicks, “Designing a mirror to realize a given projection,” J. Opt. Soc. Am. A **22**(2), 323–330 (2005). [CrossRef]

**35. **J. Hermann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. **70**(1), 28–35 (1980). [CrossRef]

**36. **R. J. Koshel, “Simplex optimization method for illumination design,” Opt. Lett. **30**(6), 649–651 (2005). [CrossRef]

**37. **S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by Simulated Annealing,” Science **220**(4598), 671–680 (1983). [CrossRef]

**38. **https://ww2.mathworks.cn/help/gads/how-simulated-annealing-works.html.

**39. **R. J. Koshel, “Aspects of illumination system optimization,” Proc. SPIE **5529**, 206–217 (2004). [CrossRef]

**40. **A. Hirst, J. Muschaweck, and P. Benítez, “Fast, deterministic computation of irradiance values using a single extended source in 3D,” Opt. Express **26**(14), A651–A656 (2018). [CrossRef]