## Abstract

We propose a method for designing optical elements with two freeform refracting surfaces generating prescribed non-axisymmetric irradiance distributions in the case of an extended light source. The method is based on the representation of the optical surfaces as bicubic splines and on the subsequent optimization of their parameters using a quasi-Newton method. For the fast calculation of the merit function, we propose an efficient version of the ray tracing method. Using the proposed approach, we design optical elements generating uniform square-shaped irradiance distributions in the far- and near-field. The designed elements are very compact (the height-to-source ratio is only 1.6) and, while providing a high lighting efficiency of 89%, generate highly uniform distributions (the ratio between minimum and average irradiance values in the prescribed square-shaped region exceeds 0.9).

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

## Corrections

3 December 2020: A typographical correction was made to the acknowledgments section.

## 1. Introduction

Nowadays, light emitting diodes (LEDs) have great advantages in practical lighting applications owing to their long lifetime, high efficiency and compact size as compared to other light sources. Secondary optics such as freeform refractive optical elements is often employed to redirect the light emitted by the LED source and to generate a prescribed light pattern. In general, the problem of calculating a refractive optical element generating a prescribed light distribution belongs to the class of inverse problems of nonimaging optics. For the case of the point-source approximation, in which the LED die is considered as a point light source, a large number of efficient design methods have been proposed [1–10]. At the same time, the design of compact secondary optics operating with high-power extended LEDs is currently of great interest. In such design problems, the element-to-LED-size ratio can be rather small (of about two or even less), so that the spatial extent of the LED source cannot be disregarded. In this case, the point-source solution is usually very different from the optimal one.

Inverse problems of designing optical elements operating with extended light sources are significantly more complex from the theoretical point of view than the corresponding point-source problems. For solving these problems, only a few approaches have been proposed, such as the feedback design techniques [11–13], “direct optimization” methods [14–17] and the methods that can be referred to as “edge-ray methods” (since they utilize the edge-ray principle for controlling the generated distribution) [18–24]. The feedback design methods use a point-source solution as the initial design, which is then iteratively modified to adapt to a real extended light source. In direct optimization methods, the surface of the optical element is represented in a certain parametric form, e.g. as a polynomial [14] or a spline [15,17]. The parameters defining the surface are considered as arguments of a merit function describing the deviation of the generated distribution from the required one. Then, the calculation of the optical surface is reduced to the minimization of the merit function with the use of a certain optimization algorithm. Let us note that the minimization of the merit function requires multiple calculations of the irradiance distribution generated by the element at a current set of parameters. The generated irradiance distribution is usually calculated using the ray tracing method.

The methods utilizing the edge-ray principle are based on the design techniques initially proposed for solar concentrators [18–20]. In recent years, several efficient versions of the edge-ray methods for designing axisymmetric double-surface refractive optical elements were developed [21–25]. In these methods, the profiles of the optical surfaces are constructed sequentially (segment by segment) using the edge rays in such a way so that a required relation between the light flux in the current fragment of the target distribution and the “lit area” of the optical element illuminating this fragment is provided. In a paper recently published in Optica, the so-called wavefront tailoring method (WTM) for designing optical elements generating non-axisymmetric far-field irradiance distributions was presented [26]. The WTM is a further development of the SMS3D design methods [1,27,28] and is based on the transformation of the “input” wavefronts from the edge points of the light source to the “output” wavefronts providing the generation of the required distribution. An example of calculating a highly compact double-surface optical element generating a nearly perfect square-shaped irradiance distribution presented in Ref. [26] demonstrates high performance of the WTM. Therefore, Ref. [26] can be considered as the most significant result in the field of designing optical elements operating with an extended light source. At the same time, in the opinion of the present authors, the WTM is a quite sophisticated method in terms of its understanding, software implementation, and practical use. In addition, the WTM has a number of significant limitations: in particular, it is unclear, whether the method can be extended to the cases of a non-square illuminated region and a non-Lambertian light source.

In the present work, we propose a direct optimization method for solving inverse problems of designing double-surface refractive optical elements. Direct optimization methods are often criticized for high computational complexity caused by the need to solve the direct problem (i.e. the problem of calculating the irradiance distribution generated by the element) multiple times [1, 21,22]. An efficient version of the ray-tracing method presented in this work enables solving the direct problem in fractions of a second on a common desktop PC. As it is shown below, this, combined with the use of a quasi-Newton optimization method implemented in MATLAB, allowed the authors to design an optical element with the characteristics exceeding the ones of the element presented in Ref. [26]. Moreover, the proposed method, as compared to the WTM, is more versatile and does not suffer from the limitations imposed by the far field condition. As an example demonstrating it, we present a highly compact optical element generating a square-shaped irradiance distribution in the near field and demonstrate the possibility of using this element for a thin direct-lit backlight system.

## 2. Problem statement

Let us consider the problem of calculating a double-surface refractive optical element. Assume that in a domain *S* of the plane $z = 0$, an extended light source radiating into the hemisphere $z > 0$ is located. The center of the source (of the domain $S$) coincides with the origin of coordinates. The radiation emitted by the source is described by the brightness function $L({{\mathrm{\varphi}} ,{\mathrm{\psi}} ,{\mathbf x}} ),\,{\mathbf x} \in S$, where ${\mathrm{\varphi}} \in [{0,2\pi } )$ and ${\mathrm{\psi}} \in [{0,{\pi / 2}} )$ are the polar and azimuthal angles of the spherical coordinate system, respectively, and ${\mathbf x} = ({x,y} )$ are the Cartesian coordinates in the source plane. The emitted radiation impinges on an optical element with a refractive index $n > 1$ having two refracting surfaces (Fig. 1). The optical surfaces are defined in spherical coordinates by the lengths of the radius vectors ${r_m}({{\mathrm{\varphi}} ,{\mathrm{\psi}} } ),\,\,m = 1,2$. The refractive index of the surrounding medium over and under the element equals ${n_0} = 1$. The problem consists in calculating the surfaces of the element ${r_m}({{\mathrm{\varphi}} ,{\mathrm{\psi}} } ),\,\,m = 1,2$ from the condition of generating a prescribed irradiance distribution ${E_t}({\mathbf u} ),\,{\mathbf u} \in D\,$ in the plane $z = f$, where ${\mathbf u} = ({u,v} )$ are the Cartesian coordinates in this plane.

## 3. Solution of the direct problem

As it was mentioned above, optimization methods for solving the inverse problem of designing an optical element generating a prescribed irradiance distribution include the multiple solution of the direct problem, which consists in calculating the irradiance distribution $E({\mathbf u} )$ generated by an optical element with certain surfaces ${r_m}({{\mathrm{\varphi}} ,{\mathrm{\psi}} } ),\,\,m = 1,2$. In this section, we discuss the method for solving the direct problem used in the present work. To calculate $E({\mathbf u} )$, let us first consider the calculation of the irradiance distribution $E({{\mathbf u};{{\mathbf x}_0}} )$ generated by an element $\textrm{d}S$ of the extended light source with the element center at a point ${{\mathbf x}_0} \in S$. We will describe the directions of the rays emitted from the point ${{\mathbf x}_0}$ by spherical coordinates $({{\mathrm{\varphi}} ,{\mathrm{\psi}} } )$. We assume that after being refracted by the surfaces of the optical element, the rays come to the target plane $z = f$. In this case, we can say that the optical element defines a certain ray mapping ${\mathbf u} = {\mathbf u}({{\mathrm{\varphi}} ,{\mathrm{\psi}} ;{{\mathbf x}_0}} )$ connecting the coordinates $({{\mathrm{\varphi}} ,{\mathrm{\psi}} } )$ of the emitted rays with the coordinates of their intersection with the target plane (Fig. 1). In this mapping, the coordinate of a point of the source ${{\mathbf x}_0}$ is considered as a parameter. The mapping ${\mathbf u}({{\mathrm{\varphi}} ,{\mathrm{\psi}} ;{{\mathbf x}_0}} )$ defines the irradiance distribution $E({{\mathbf u};{{\mathbf x}_0}} )$ generated in the plane $z = f$ according to the light flux conservation law:

In order to calculate the “total” irradiance distribution generated by the whole extended light source, one has to integrate Eq. (3) over the source domain:

In the present work, we calculate the irradiance distribution $E({\mathbf u} )$ of Eq. (5) using the Monte Carlo method [29]. In this method, *N* realizations $\boldsymbol{\mathrm{\xi}}_k = ({{x_k},{y_k},{{\mathrm{\varphi}}_k},{{\mathrm{\psi}}_k}} ),\,k = 1,\ldots ,N$ of a random vector with independent components are generated. We consider each realization $\boldsymbol{\mathrm{\xi}}_k$ as a ray emitted from a random point ${{\mathbf x}_k} = ({{x_k},{y_k}} )\in S$ and having a random direction $\boldsymbol{\mathrm{\alpha}}_k = ({{{\mathrm{\varphi}}_k},{{\mathrm{\psi}}_k}} )\in \Theta $. In the general case, the probability density function of the random vector $\boldsymbol{\mathrm{\xi}} = ({\mathbf x},\boldsymbol{\mathrm{\alpha}} )$ has to be chosen in a form matched to the brightness function of the source. However, for the sake of simplicity, let us assume that the coordinates of the points ${{\mathbf x}_k}$ are the realizations of a random vector uniformly distributed over the domain *S*. The directions of the rays $\boldsymbol{\mathrm{\alpha}}_k$ are assumed to be the realizations of a random vector with the probability density function

*S*. One can easily see that Eq. (8) represents the sum of “elementary” light flux values $f({{{\mathbf x}_k},{{\mathrm{\varphi}}_k},{{\mathrm{\psi}}_k}} )\Delta V$ emitted by randomly chosen elements of the source in random directions. Note that these light flux values are distributed over the vicinity of the points ${{\mathbf u}_k} = {\mathbf u}({{{\mathrm{\varphi}}_k},{{\mathrm{\psi}}_k};{{\mathbf x}_k}} )$ with a Gaussian weight (4). It is known that the expectation of the random quantity ${E_N}({\mathbf u} )$ is equal to the exact value of the integral (5), and its dispersion is inversely proportional to the square root of the number of rays

*N*(i.e. tends to zero with an increase in

*N*) [29].

In order to calculate the irradiance using Eq. (8), it is necessary to trace the rays $\boldsymbol{\mathrm{\xi}}_k$. This ray tracing is required for finding the coordinates ${{\mathbf u}_k} = {\mathbf u}({{{\mathrm{\varphi}}_k},{{\mathrm{\psi}}_k};{{\mathbf x}_k}} )$ of the points of intersection of the emitted rays with the target plane. The ray tracing procedure includes finding the points of intersection of the incident rays with the optical surfaces, calculating the directions of the refracted rays and the transmission coefficient [function $T({{\mathrm{\varphi}} ,{\mathrm{\psi}} ;{\mathbf x}} )$ in Eqs. (1) and (6)]. In this work, to calculate the points of intersection of the rays with the surfaces of the optical element, the surfaces ${r_m}({{\mathrm{\varphi}} ,{\mathrm{\psi}} } ),\,\,m = 1,2$ were triangulated. Since these surfaces are defined in a parametric form, the triangulation is simply achieved by triangulating the domain of definition $\Theta $. The number and size of the triangles in the domain $\Theta $ were chosen from the condition of representing the functions ${r_m}({{\mathrm{\varphi}} ,{\mathrm{\psi}} } ),\,\,m = 1,2$ with a required accuracy. The simplest algorithm for finding the intersection of a ray with a triangulated surface consists in sequentially checking the condition of the intersection of the ray with each triangle (primitive). The computational complexity of this algorithm is $O({{N_{tr}}} )$, where ${N_{tr}}$ is the number of primitives. In order to decrease the computational complexity, the primitives are ordered in space. This approach enables efficiently excluding from search the primitives, which the incident ray certainly does not cross. In the present work, we used the so-called k-d tree data structure [30] for binary space-partitioning, which makes it possible to use the binary search algorithm [31,32] for finding the primitive intersected by the ray with computational complexity $O({\log {N_{tr}}} )$. The direction of the refracted ray is then found using the Snell’s law. For calculating $T({{\mathrm{\varphi}} ,{\mathrm{\psi}} ;{\mathbf x}} )$, which is equal to the product of the transmission coefficients of two surfaces of the element, Fresnel equations for unpolarized light were used.

Thus, the method for solving the direct problem includes tracing *N* randomly generated rays $\boldsymbol{\mathrm{\xi}}_k = ({{\mathbf x}_k},\boldsymbol{\mathrm{\alpha}}_k)$ in order to find the coordinates ${{\mathbf u}_k} = {\mathbf u}({{{\mathrm{\varphi}}_k},{{\mathrm{\psi}}_k};{{\mathbf x}_k}} )$, and, then, calculating the generated irradiance distribution using Eq. (8).

## 4. Solution of the inverse problem

Let us now move to the problem of calculating the optical surfaces ${r_m}({{\mathrm{\varphi}} ,{\mathrm{\psi}} } ),\,\,m = 1,2$ from the condition of generating a prescribed irradiance distribution ${E_t}({\mathbf u} ),\,{\mathbf u} \in D\,$. In this work, we propose to define the functions ${r_m}({{\mathrm{\varphi}} ,{\mathrm{\psi}} } ),\,\,m = 1,2$ as bicubic splines [33] defined on the grids ${\Lambda _m} = ({{{\mathrm{\varphi}}_{m,i}},{{\mathrm{\psi}}_{m,j}}} ),i = 1,\ldots ,{N_m},j = 1,\ldots ,{M_m}$. In each grid cell ${C_{m,ij}} = \{{ {({{\mathrm{\varphi}} ,{\mathrm{\psi}} } )} |{{\mathrm{\varphi}}_{m,i}} \le {\mathrm{\varphi}} \le {{\mathrm{\varphi}}_{m,i + 1}},\,\,{{\mathrm{\psi}}_{m,j}} \le {\mathrm{\psi}} \le {{\mathrm{\psi}}_{m,j + 1}}} \}$, such a spline is a polynomial of the form

Let us denote by ${{\mathbf c}_m} = ({{r_{m,\,i,j}},\,r_{m,\,i,j}^{({\mathrm{\varphi}} )},\,\,r_{m,\,i,j}^{({\mathrm{\psi}} )},\,r_{m,\,i,j}^{({{\mathrm{\varphi}} {\mathrm{\psi}} } )}} ),\,\,m = 1,2$ the vectors consisting of the spline parameters describing the first and the second surfaces. The irradiance distribution generated by the optical element in the target plane $z = f$ depends on the spline parameters ${\mathbf c} = ({{\mathbf c}_1},{{\mathbf c}_2})$. Therefore, we will denote this distribution by $E({{\mathbf u};{\mathbf c}} )$.

If the optical surfaces are defined as bicubic splines, the problem of calculating an optical element can be formulated as a problem of minimizing a certain function $\varepsilon ({\mathbf c} )$ of the spline parameters ${\mathbf c} = ({{{\mathbf c}_1},{{\mathbf c}_2}} )$, which characterizes the deviation of the irradiance distribution generated by the element $E({{\mathbf u};{\mathbf c}} )$ from the required distribution ${E_t}({\mathbf u} ),\,{\mathbf u} \in D\,$. Here, we use the root-mean-square deviation (RMSD) of the generated distribution from the required one. In this case, the design of the optical element is reduced to the following optimization problem:

*fminunc*MATLAB routine. In the simplest case, in order to use this routine, it is sufficient to calculate only the function (10) (in this case, the estimates of the gradient and Hessian of this function are calculated numerically using finite-difference expressions). The calculation of $\varepsilon ({\mathbf c} )$ is reduced to calculating the irradiance distribution $E({{\mathbf u};{\mathbf c}} )$ generated by the element with the given parameters vector ${\mathbf c}$. The method of calculating the distribution $E({{\mathbf u};{\mathbf c}} )$ is described in the previous section and was also implemented in MATLAB.

Thus, the optical element is calculated using the following iterative optimization procedure:

- 1. The parameters of the design problem (size and the brightness distribution of the light source, required irradiance distribution, material of the optical element, distance to the target plane, etc.) are specified.
- 2. The initial optical surfaces ${r_m}({{\mathrm{\varphi}} ,{\mathrm{\psi}} } ),\,\,m = 1,2$ are defined and approximated by bicubic splines. The initial spline parameters ${\mathbf c} = ({{\mathbf c}_1},{{\mathbf c}_2})$ are calculated.
- 3. The irradiance distribution $E({{\mathbf u};{\mathbf c}} )$ generated by the element with the given parameters vector ${\mathbf c}$ is calculated using the ray-tracing procedure described in Section 4.
- 4. The RMSD $\varepsilon ({\mathbf c} )$ is calculated. If $\varepsilon ({\mathbf c} )< {\varepsilon _{req}}$, where ${\varepsilon _{req}}$ is the required accuracy of the target irradiance distribution, the current solution is used and the iterative process is stopped. Otherwise, continue to the next step.
- 5. The vector of spline parameters is corrected using quasi-Newton BFGS method implemented in the
*fminunc*MATLAB routine. Return to step 3.

Let us note that for representing the surfaces of the optical element, instead of bicubic splines one could choose other functions such as polynomials [14], B-splines [9,35] or nonrational B-splines [14,36]. The authors decided to use bicubic splines (9) defined in spherical coordinates, since this representation proved efficient for solving various inverse problems of designing optical elements for a point light source [17,37,38].

## 5. Design examples

#### 5.1 Double-surface element generating a far-field irradiance distribution

As a first example, we designed an optical element (refractive index $n = 1.491$) generating a “smoothed uniform” square-shaped irradiance distribution in the plane $z = f = 1000\textrm{ mm}$ for an extended Lambertian light source with the size of $1 \times 1\textrm{ m}{\textrm{m}^2}$. The required distribution is constant in a square-shaped region with the side $d = 2f \cdot tg({30^\circ } )\approx 1155\,\textrm{mm}$ and then gradually decreases to zero [ Fig. 2(c)]. The cross-section of the target distribution is shown with a dashed line in Fig. 3(d). Such a distribution was chosen for comparison with the results of the recent paper [26]. In that paper, an optical element is designed using the wavefront tailoring method, which generates a similar irradiance distribution for an extended source with the same size. In what follows, we will refer to this element as WTM-element. The WTM-element is very compact: the ratio of the element height to source size (diagonal of the square-shaped LED die) is equal to 1.67. This is a very good result, taking into account the high quality of the irradiance distribution generated by this element [26].

Let us discuss in detail the design of an element generating the described irradiance distribution using the proposed method presented above. Due to the symmetry of the generated irradiance distribution, the designed element possesses the following planes of symmetry: $y = 0$ $({\mathrm{\varphi}} = 0)$, $x = 0$ $({\mathrm{\varphi}} = {\pi / 2})$, and $y ={\pm} x$ (${\mathrm{\varphi}} = {\pi / 4}$ and ${\mathrm{\varphi}} = {{3\pi } / 4}$). Therefore, it is sufficient to calculate the optical surfaces only at ${\mathrm{\varphi}} \in [{0,\,{\pi / 4}} ]$, whereas the “complete” surfaces of the optical element (at ${\mathrm{\varphi}} \in [{0,\,2\pi } )$) can be obtained from the surface fragments calculated at ${\mathrm{\varphi}} \in [{0,\,{\pi / 4}} ]$ by sequential reflection with respect to the symmetry planes ${\mathrm{\varphi}} = {\pi / 4}$, ${\mathrm{\varphi}} = {\pi / 2}$, and ${\mathrm{\varphi}} = 0$. Thus, as optimization parameters (parameters of bicubic splines defining the radius vectors of the element surfaces), it suffices to consider the spline parameters for ${\mathrm{\varphi}} \in [{0,\,{\pi / 4}} ]$. According to Section 4, the parameters of the splines are found by solving the optimization problem (10). As a starting point of the optimization process, we used an axisymmetric element designed for a point light source and generating a uniform irradiance distribution in a circle with radius $R = 690\textrm{ mm}$ in the plane $z = 1000\textrm{ mm}$. The calculation of such elements was described in detail in Refs. [38,39] and is reduced to the integration of three explicit ordinary differential equations. The surfaces of the calculated “starting-point” element were approximated by bicubic splines. This optical element is shown in Fig. 2(а). The spline nodes to be optimized, which represent the inner and outer surfaces of the element at ${\mathrm{\varphi}} \in [{0,\,{\pi / 4}} ]$, are shown in Fig. 2(a) with blue and red dots. Note that in order for the “complete” optical surfaces to have no sharp bends (derivative discontinuities), the derivatives of the splines $\,r_m^{({\mathrm{\varphi}} )}$ and $\,r_m^{({{\mathrm{\varphi}} {\mathrm{\psi}} } )}$ at the spline grid nodes corresponding to ${\mathrm{\varphi}} = 0$ and ${\mathrm{\varphi}} = {\pi / 4}$ were set equal to zero and excluded from the optimization process. At the nodes located at the vertices of the first and second surfaces (at ${\mathrm{\psi}} = 0$), both first- and second-order derivatives were set equal to zero, and only the lengths of the radius vectors ${r_m}$ were considered as optimization parameters. Under these conditions and for the chosen spline grids, the number of the optimized spline parameters amounts to 13 and 85 for the inner and outer surfaces, respectively. Figure 2(b) shows the irradiance distribution generated by the optical element in Fig. 2(a) when illuminated by an extended Lambertian light source with the dimensions of $1 \times 1\textrm{ m}{\textrm{m}^2}$. This distribution is very different from the required distribution shown in Fig. 2(c), however, as we will show below, the optimization of the chosen initial surfaces enables obtaining a good solution to the problem.

The optimization of the optical element [the minimization of the function (10)] was performed using a standard MATLAB routine *fminunc*. For the calculation of the function $E({{\mathbf u};{\mathbf c}} )$ during the optimization, ray tracing technique described in Section 3 was used. In each calculation, $N = 500000$ rays were traced. This number of rays provides an acceptable tradeoff between the computation time and accuracy (with a twofold increase in the number of rays, the obtained irradiance distribution $E({{\mathbf u};{\mathbf c}} )$ changes by no more than 2–3%). At $N = 500000$, the average computation time of the function $\varepsilon ({\mathbf c} )$ on a desktop PC (Intel Core i7-3770) amounted to 0.67 seconds. During the optimization process, the function $\varepsilon ({\mathbf c} )$ was evaluated 19296 times, and the total calculation time was of about 3.5 hours. The resulting optical element is shown in Fig. 3(a). The height-to-source ratio for the designed element [see Fig. 3(b)] is equal to 1.6, which is slightly better than for the WTM-element [26], which generates a similar irradiance distribution. Besides, the width of the designed element is also less than the width of the WTM-element (compare the black cross-section in Fig. 3(b) with the upper panel of Fig. 5(a) in Ref. [26]). It is important to note that in contrast to the element designed in the present work [Figs. 3(a) and 3(b)], the WTM-element is less suited for the fabrication and practical use. Indeed, that element consists of two separate surface fragments that do not fully encompass the light source and therefore “capture” only a part of the light flux emitted by the source (see Fig. 5 in Ref. [26]).

The irradiance distribution generated by the optimized element is shown in Fig. 3(c). Let us note that in order to validate our calculations, the distribution in Fig. 3(c) was calculated using the commercial ray-tracing software TracePro [40]. The obtained distribution is in good agreement with the required one, the cross-section of which is shown with a dashed line in Fig. 3(d). Following Ref. [26], in order to characterize the uniformity of the generated distribution, we use the ratio ${{{E_{\min }}} / {{E_{average}}}}$ between the minimum and average irradiance in the central “flat-top” region. For the distribution shown in Fig. 3(c), ${{{E_{\min }}} / {{E_{average}} = 0.92}}$. For the WTM-element generating a similar distribution, the value ${{{E_{\min }}} / {{E_{average}}}}$ is slightly less and amounts to 0.9 [26]. In addition, the optical element calculated in this work has a higher lighting efficiency, which is defined as a fraction of the light flux emitted by the source reaching the target plane. For the WTM-element, the lighting efficiency, which was calculated without taking into account the Fresnel losses, was 80% [26]. For the optical element shown in Fig. 3, the lighting efficiency calculated using TracePro taking into account these losses, amounts to 89.1%. The achieved efficiency is quite close to the theoretical limit. Indeed, for a plane interface between air and a material with refractive index *n*, the Fresnel transmission coefficient at normal incidence is equal to ${T_1} = {{4n} / {{{({n + 1} )}^2}}}$. Since the optical element has two refractive surfaces, the resulting transmission coefficient can be estimated as $T \le T_1^2 = {[{{{4n} / {{{({n + 1} )}^2}}}} ]^2}$. In this case, the upper estimate of the theoretically achievable lighting efficiency is $\eta = 1 - {T^2}$. In the considered example, $n = 1.491$ and this upper estimate amounts to 92.4%.

#### 5.2 Double-surface element generating a near-field irradiance distribution

The example presented above demonstrates high performance of the proposed optimization method when generating a uniform irradiance distribution in the far field, when the dimensions of the optical element can be neglected in comparison with the distance to the target plane. It is important to note that the proposed method, as compared to the method of Ref. [26], is more versatile and can also be used for designing optical elements generating prescribed irradiance distributions in the near field. Such problems are particularly relevant in the design of LED backlight systems.

To illustrate the capabilities of the method in the near-field case, we designed an optical element generating a smoothed uniform square-shaped irradiance distribution in the near plane $z = f = 10\textrm{ mm}$ in the case of an extended Lambertian light source with the dimensions of $1 \times 1\textrm{ m}{\textrm{m}^2}$. Similarly to the previous example, the target irradiance is constant in a square region with the side $d = 50\textrm{ mm}$ and then gradually decreases to zero [Fig. 4(c)]. The cross-section of the target distribution is shown with a dashed line in Fig. 5(d).

The calculation of the “near-field” optical element was completely similar to the previous example. Again, due to the symmetry of the generated distribution, we optimized the parameters of the bicubic splines representing the optical surfaces at ${\mathrm{\varphi}} \in [{0,\,{\pi / 4}} ]$ only. As a starting point of the optimization process, we used an axisymmetric element designed for a point light source and generating a uniform irradiance distribution in a circle with radius $R = 24\textrm{ mm}$ in the near plane $z = 10\textrm{ mm}$ [38]. The surfaces of the calculated “starting-point” near-field element are shown in Fig. 4(а). The spline nodes to be optimized, which represent the inner and outer surfaces of the element at ${\mathrm{\varphi}} \in [{0,\,{\pi / 4}} ]$, are shown in Fig. 4(a) with blue and red dots. Figure 4(b) shows the irradiance distribution generated by the optical element in Fig. 4(a) in the case of an extended Lambertian light source with the dimensions of $1 \times 1\textrm{ m}{\textrm{m}^2}$. As in the previous example, this distribution is very different from the required distribution shown in Fig. 4(c).

During the optimization, the function $\varepsilon ({\mathbf c} )$ was calculated 7250 times, and the total calculation time was of about 2 hours. The obtained near-field optical element is shown in Fig. 5(a). As in the previous example, the height-to-source ratio for the calculated element [see Fig. 5(b)] is equal to 1.6. The generated irradiance distribution is shown in Fig. 5(c). As in the first example, we calculated this distribution using the commercial ray-tracing software TracePro. The obtained distribution has a high uniformity ${{{E_{\min }}} / {{E_{average}}}} = 0.94$ in the central flat-top area and is in good agreement with the target irradiance distribution [see Fig. 4(d)]. The lighting efficiency of the designed element is high and amounts to 89.0%. Thus, in the near-field case, the proposed method also enables calculating optical elements with very high performance characteristics.

The designed near-field element can be used in the so-called direct-lit backlight systems. Indeed, using a set of optical elements located at the nodes of a square grid and directly illuminating square-shaped regions, it is possible to uniformly illuminate a given region with a large area. A uniform irradiance distribution in this case can be achieved by overlapping of the irradiance distributions generated by the adjacent optical elements. To illustrate this possibility, we calculated using TracePro the irradiance distribution generated by a set of five optical elements arranged in the form of a cross. In the simulation, the light sources were located at the following five points of the plane $z = 0$: $({0,0} )$, $({ \pm 60\textrm{ mm},0} )$, and $({0, \pm 60\textrm{ mm}} )$. Accordingly, five identical optical elements were placed above the sources. The irradiance distribution generated by this set of elements is shown in Fig. 6 and confirms the possibility of uniformly illuminating a larger area.

In order to characterize the compactness of a backlight system, the distance-to-height ratio (DHR) is often used, defined as the ratio of the distance between two adjacent light sources (LEDs) to the thickness of the backlight system (in the considered case, the distance to the target plane $z = 10\textrm{ mm}$). In the presented example, the DHR value is equal to 6, which, to the best of our knowledge, is a record value for modern direct-lit backlight systems [41–44].

## 6. Conclusion

We presented an optimization method for designing optical elements with two freeform surfaces operating with an extended light source and generating prescribed irradiance distributions. In the method, the surfaces of the optical element are represented as bicubic splines. The spline parameters are considered as arguments of a merit function describing the root-mean-square deviation of the generated distribution from the required one. For fast calculation of the merit function, we proposed an efficient version of the ray tracing method. The minimization of the merit function was performed using a quasi-Newton method implemented in MATLAB.

Using the proposed method, we designed compact optical elements generating smoothed uniform square-shaped irradiance distributions in the far- and near-field in the case of an extended Lambertian light source. The ratio between the height of the designed elements and the size of the light source was only 1.6. Despite such a compact size, the lighting efficiency of the calculated elements amounts to 89%, and the uniformity of the generated distributions, defined as the ratio between the minimum and average irradiance, exceeds 0.9. To the best of our knowledge, the performance of the designed elements is currently record-breaking.

In the present work, we limit our consideration to the design of optical elements generating only square-shaped irradiance distributions. However, we believe that the proposed optimization method can also handle irradiance distributions having more complex shapes. Moreover, for representing the surfaces of the optical element, instead of bicubic splines, one could use other functions such as B-splines or nonrational B-splines. The study of the optimization method performance for more complex irradiance distributions and for different surface representations will be a topic of a further research.

## Funding

Russian Science Foundation (18-19-00326); Russian Foundation for Basic Research (18-07-00982, 18-29-03067); Ministry of Science and Higher Education of the Russian Federation (State assignment to FSRC “Crystallography and Photonics” RAS).

## Acknowledgment

This work was funded by Russian Science Foundation (the method of solving the inverse problem in Section 4 and the design of the optical elements in Section 5), by Russian Foundation for Basic Research (the method of solving the direct problem in Section 3), and by the Ministry of Science and Higher Education of Russian Federation (implementation of the surface triangulation and k-d tree data structure in the ray tracing simulation software).

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **R. Wu, Z. Feng, Z. Zheng, R. Liang, P. Benítez, J. C. Miñano, and F. Duerr, “Design of freeform illumination optics,” Laser Photonics Rev. **12**(7), 1700310 (2018). [CrossRef]

**2. **R. Wu, L. Xu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “Freeform illumination design: a nonlinear boundary problem for the elliptic Monge–Ampére equation,” Opt. Lett. **38**(2), 229–231 (2013). [CrossRef]

**3. **X. Mao, S. Xu, X. Hu, and Y. Xie, “Design of a smooth freeform illumination system for a point light source based on polar-type optimal transport mapping,” Appl. Opt. **56**(22), 6324–6331 (2017). [CrossRef]

**4. **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]

**5. **K. Brix, Y. Hafizogullari, and A. Platen, “Designing illumination lenses and mirrors by the numerical solution of Monge–Ampère equations,” J. Opt. Soc. Am. A **32**(11), 2227–2236 (2015). [CrossRef]

**6. **N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Computation of double freeform optical surfaces using a Monge–Ampère solver: Application to beam shaping,” Opt. Commun. **439**, 251–259 (2019). [CrossRef]

**7. **Y. Schwartzburg, R. Testuz, A. Tagliasacchi, and M. Pauly, “High-contrast computational caustic design,” ACM Trans. Graph. **33**(4), 1–11 (2014). [CrossRef]

**8. **V. Oliker, “Controlling light with freeform multifocal lens designed with supporting quadric method (SQM),” Opt. Express **25**(4), A58–A72 (2017). [CrossRef]

**9. **L. L. Doskolovich, D. A. Bykov, E. S. Andreev, E. A. Bezus, and V. Oliker, “Designing double freeform surfaces for collimated beam shaping with optimal mass transportation and linear assignment problems,” Opt. Express **26**(19), 24602–24613 (2018). [CrossRef]

**10. **L. L. Doskolovich, D. A. Bykov, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions,” Opt. Express **27**(9), 13083–13097 (2019). [CrossRef]

**11. **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]

**12. **Z. Li, S. Yu, L. Lin, Y. Tang, X. Ding, W. Yuan, and B. Yu, “Energy feedback freeform lenses for uniform illumination of extended light source LEDs,” Appl. Opt. **55**(36), 10375–10381 (2016). [CrossRef]

**13. **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]

**14. **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]

**15. **F. Fournier and J. Rolland, “Optimization of freeform lightpipes for light-emitting-diode projectors,” Appl. Opt. **47**(7), 957–966 (2008). [CrossRef]

**16. **S. Zhao, K. Wang, F. Chen, Z. Qin, and S. Liu, “Integral freeform illumination lens design of LED based pico-projector,” Appl. Opt. **52**(13), 2985–2993 (2013). [CrossRef]

**17. **M. A. Moiseev and L. L. Doskolovich, “Design of refractive spline surface for generating required irradiance distribution with large angular dimension,” J. Mod. Opt. **57**(7), 536–544 (2010). [CrossRef]

**18. **H. Ries and R. Winston, “Tailored edge-ray reflectors for illumination,” J. Opt. Soc. Am. A **11**(4), 1260–1264 (1994). [CrossRef]

**19. **W. T. Welford and R. Winston, * High Collection Nonimaging Optics* (Academic Press, 1989).

**20. **J. Chaves, * Introduction to Nonimaging Optics*, 2nd ed. (CRC Press, 2016).

**21. **R. Wu, C. Y. Huang, X. Zhu, H.-N. Cheng, and R. Liang, “Direct three-dimensional design of compact and ultra-efficient freeform lenses for extended light sources,” Optica **3**(8), 840–843 (2016). [CrossRef]

**22. **R. Wu, H. Hua, P. Benítez, and J. C. Miñano, “Direct design of aspherical lenses for extended non-Lambertian sources in two-dimensional geometry,” Opt. Lett. **40**(13), 3037–3040 (2015). [CrossRef]

**23. **S. Hu, K. Du, T. Mei, L. Wan, and N. Zhu, “Ultra-compact LED lens with double freeform surfaces for uniform illumination,” Opt. Express **23**(16), 20350–20355 (2015). [CrossRef]

**24. **X. Li, P. Ge, and H. Wang, “Prescribed intensity in 3D rotational geometry for extended sources by using a conversion function in 2D design,” Appl. Opt. **56**(6), 1795–1798 (2017). [CrossRef]

**25. **X. Li, P. Ge, and H. Wang, “An efficient design method for LED surface sources in three-dimensional rotational geometry using projected angle difference,” Lighting Res. Technol. **51**(3), 457–464 (2019). [CrossRef]

**26. **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]

**27. **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]

**28. **O. Dross, R. Mohedano, P. Benítez, J. C. Miñano, J. Chaves, J. Blen, M. Hernández, and F. Muñoz, “Review of SMS design methods and real-world applications,” Proc. SPIE **5529**, 35 (2004). [CrossRef]

**29. **C.P. Robert and G. Casella, * Monte Carlo Statistical Methods*, 2nd ed. (Springer, 2004).

**30. **J. L. Bentley, “Multidimensional binary search trees used for associative searching,” Commun. ACM **18**(9), 509–517 (1975). [CrossRef]

**31. **T. H. Cormen, C. E. Leiserson, and R. L. Rivest, * Introduction to Algorithms* (MIT Press, 1989).

**32. **R. Sedgewick, * Algorithms in C++* (Addison-Wesley, 1992).

**33. **C. De Boor, * A Practical Guide to Splines* (Springer, 2001).

**34. **Bicubic interpolation, https://en.wikipedia.org/wiki/Bicubic_interpolation.

**35. **R. Wu, Z. Zheng, H. Li, and X. Liu, “Constructing optical freeform surfaces using unit tangent vectors of feature data points,” J. Opt. Soc. Am. A **28**(9), 1880–1888 (2011). [CrossRef]

**36. **R. Wu, J. Sasián, and R. Liang, “Algorithm for designing free-form imaging optics with nonrational B-spline surfaces,” Appl. Opt. **56**(9), 2517–2522 (2017). [CrossRef]

**37. **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]

**38. **S. V. Kravchenko, M. A. Moiseev, E. V. Byzov, and L. L. Doskolovich, “Design of axisymmetric double-surface refractive optical elements generating required illuminance distributions,” Opt. Commun. **459**, 124976 (2020). [CrossRef]

**39. **S. V. Kravchenko, M. A. Moiseev, L. L. Doskolovich, and N. L. Kazanskiy, “Design of axis-symmetrical optical element with two aspherical surfaces for generation of prescribed irradiance distribution,” Computer Opt. **35**, 467–472 (2011). [CrossRef]

**40. **TracePro, “Illumination and non-imaging optical design & analysis tool,” https://www.lambdares.com/tracepro.

**41. **C. Leiner, W. Nemitz, F. P. Wenzl, and C. Sommer, “Ultrathin free-form micro-optical elements for direct-lit applications with a large distance-height ratio,” OSA Continuum **1**(4), 1144–1157 (2018). [CrossRef]

**42. **R. Wu, Z. Zheng, H. Li, and X. Liu, “Optimization design of irradiance array for LED uniform rectangular illumination,” Appl. Opt. **51**(13), 2257–2263 (2012). [CrossRef]

**43. **J. Zheng and K. Qian, “Designing single LED illumination distribution for direct-type backlight,” Appl. Opt. **52**(28), 7022–7027 (2013). [CrossRef]

**44. **L. L. Doskolovich, O. I. Petrova, and M. A. Moiseev, “LED lighting system based on modules that form uniform irradiance of a hexagonal region,” J. Opt. Technol. **78**(2), 105–109 (2011). [CrossRef]