Initial design with L<sup>2</sup> Monge-Kantorovich theory for the Monge–Ampère equation method in freeform surface illumination design

Rengmao Wu; Yaqin Zhang; Mohamed M. Sulman; Zhenrong Zheng; Pablo Benítez; Juan C. Miñano

doi:10.1364/OE.22.016161

1. Introduction

Illumination design, which controls the distribution of light from a source to produce a target illumination, is a classical problem of nonimaging optics [1]. The traditional method of solving this problem is using the quadric surfaces, such as paraboloid, ellipsoid and spherical surface, to direct the light. The quadric surfaces, though, have been playing an important role in illumination design, the low degree of design freedom of the quadric surfaces, is a frustrating problem faced by the designer. With the increasing emphasis on energy conservation and practical issues, it is becoming a preferred route to solve the problem of illumination design with freeform surface which has higher degree of design freedom [2–11]. The difficulty in the freeform surface illumination design is how to exactly tailor the freeform surface only on the basis of the given intensity distribution of the source and the target illumination. A variety of methods have been developed to solve this problem in the past two decades. The Simultaneous Multiple Surfaces method proposed by Benítez and Miñano allows to simultaneously construct two freeform surfaces that exactly transform two input wavefronts into two output wavefronts [2]. For complex designs, the difficulty of this method lies in establishing a specific relationship between the input and output wavefronts. Optimization method is an iteration of successively finding some appropriate variable values to reduce the merit function value with certain optimization algorithms [3,4]. Due to the Monte-Carlo raytracing and the limited number of optimization variables, the optimization method may not be the best choice in illumination design. The first-order Partial Differential Equation (PDE) method is popular in solving this problem due to its ease of use and fast computation [5,6]. The continuity of the freeform surface designed by the PDE method, however, depends on the integrability of the mapping, and Fournier, using the Supporting Ellipsoids method developed by Oliker [7], has shown the difficulty in obtaining such an integrable mapping [8]. The optimal mass transport also provides an alternative way to get a mapping for freeform surface design [9,10]. Since such a mapping does not include any information about the freeform surface, the mapping may not be integrable and the limitations of these methods in illumination design still need to be further explored [9,10]. The Monge–Ampère (MA) equation method proposed by Ries has shown its elegance in solving the problem of freeform surface illumination [11]. Since the numerical solution of this method was not detailed, this MA method has been covered with a veil of mystery for a long time.

In our previous work [12], we convert the problem of freeform surface illumination design into a nonlinear boundary problem for the elliptic MA equation based on ideal source assumption and introduce an approach to find the numerical sollution of this mathematial problem. By this MA method, the boundary incident rays are mapped to the boundary of the target pattern, and the inner incident rays are pushed to satisfy the elliptical MA equation which describes the conservation and redistribution of energy during the beam-shaping process. Taking the collimated beam shaping for example, the mathematical model of the single freeform surface design is given by [13]

{\begin{cases} A_{1} (z_{x x} z_{y y} - z_{x y}^{2}) + A_{2} z_{x x} + A_{3} z_{y y} + A_{4} z_{x y} + A_{5} = 0 \\ B C : {\begin{cases} t_{x} = t_{x} (x, y, z, z_{x}, z_{y}) \\ t_{y} = t_{y} (x, y, z, z_{x}, z_{y}) \end{cases} : \partial S_{1} \to \partial S_{2} \end{cases}

where, the coefficients which can be written as A_i = A_i(x, y, z, z_x, z_y) (i = 1,…,5) are functions of x, y, z, z_x and z_y; t_x and t_y are the x- and y-coordinates of the target point on the illumination plane; ∂S₁ and ∂S₂ represent the boundaries of the domains S₁ of the incident beam and S₂ of the target pattern, respectively. Then, a numerical approach based on the Newton’s method is developed to solve this design model (Here, we assume F₁(X) = 0 represents the nonlinear equations of Eq. (1) obtained by the numerical approach.) [13]. The MA method has shown both mathematical and practical interests in freeform surface illumination design for compact light source and collimated beam [13,14]. The Newton’s method which has quadratic convergence is a method for finding better approximations to the roots of a real-valued function. However, for a highly nonlinear function, the convergence of Newton’s method is strongly determined by the initial value. The coefficients A_i (i = 1,…,4) of the MA equation given in Eq. (1), are usually highly nonlinear in z_x and z_y. Consequently, the convergence of Eq. (1) is strongly determined by the initial design, as illustrated in Fig. 1.The initial design has become one of the most urgent issues to be solved of the MA method in illumination design.

Fig. 1 Some practical considerations for the convergence of the MA method. x₀ and x₁, respectively, represent the initial value and the x-intercept obtained from the first iteration. (a) The iteration diverges at the first iteration. In this case, x₁ may deteriorate the optical performance of the surface. For example, one will uaually encounter complications associated with singularities of the freeform surface, and the method cannot be iterated (This will be demonstrated in the first design example). (b) Assume F₁(x₀) = F₁(x₁) and |F₁’(x₀)| = |F₁’(x₁)|. The iteration will go into an infinite loop in this case. (c) x₀ may be a good initial guess of a pure mathematical problem in this case. The convergence of the MA method, however, still depends on the optical performance of the freeform surface. For a freeform lens design, the MA method can converge successfully only and if only the total internal reflection will not take place on the freeform surface in the initial design (This will be demonstrated in the second design example).

Download Full Size | PDF

In order to improve the MA method and promote its application in illumination design, we address the initial design of the MA method in this paper on the basis of the L² Monge-Kantorovich (LMK) theory. A new approach of computing a numerical solution of the LMK problem is presented here, and the optimal mapping of this problem is found by solving an equivalent problem not relying on the gradient of time. The approach converges quite fast and is very flexible in imposing the boundary condition. Most importantly, the MA method is significantly improved by applying this new approach to the initial design.

The paper is organized as follows. In Section 2, we will review the LMK problem and present an efficient approach to solve the LMK problem. Then, we will introduce some characteristics of the approach and explore some limitations of the LMK theory in freeform surface illumination design in Section 3. Also in this section, comparisons will be made between the designs respectively obtained from the LMK theory and the MA method for both on-axis illumination and off-axis illumination. In Section 4, three design examples, including the beam shaping of compact light source and collimated beam, will be given to illustrate the significant improvement of the MA method achieved by the application of the LMK theory in the initial design. Finally, we will conclude with some remarks and discussions in Section 5.

2. L² Monge-Kantorovich problem and the new approach

2.1. L² Monge-Kantorovich problem

Optimal mass transport is an important problem arising in a number of applications including econometrics, image registration [15], mesh generation [16], automatic control, and reflector design [17]. However, numerical solution of this problem has proved very challenging. The optimal mass transport problem, also known as Monge-Kantorovich problem, can be defined as follows: Let ρ₀(ξ) and ρ₁(η) be two normalized non-negative density functions with compact supports Ω₀ and Ω₁, where ξ, η ∈R^d, d≥1 is the space dimension. Here the data must satisfy the condition that the total mass is conserved:

\int_{Ω_{0}} ρ_{0} (ξ) d ξ = \int_{Ω_{1}} ρ_{1} (η) d η

Assume ϕ(ξ) is a smooth one-to-one mapping that takes the density ρ₀(ξ) to the density ρ₁(η). Then, according to the coordinate transformation between ξ and η, Eq. (2) leads to

\det (\nabla ϕ (ξ)) ρ_{1} (ϕ (ξ)) = ρ_{0} (ξ)

where det(∇ϕ(ξ)) is the determinant of the Jacobian matrix of the mapping ϕ(ξ). The quantity of det(∇ϕ(ξ)) represents the expansion (or contraction) of an infinitesimal area element of Ω₀ during the mass transport. All the points contained in the infinitesimal area element will converge to a point, a straight line or a curve when det(∇ϕ(ξ)) = 0.

There are usually many such mappings which can satisfy Eq. (3), however, there is a unique optimal mapping $ϕ = \bar{ϕ}$ that minimizes the transport cost

C (ϕ) = {\int_{R^{d}} | ϕ (ξ) - ξ |}^{2} ρ_{0} (ξ) d ξ

Moreover, the optimal mapping

\bar{ϕ}

can be expressed as the gradient of a convex function w. That is,

\bar{ϕ} = \nabla w

. The infimum of Eq. (4) is the square of the L² Kantorovich-Wasserstein distance, and the optimal transport, in this case, is called the L² Monge-Kantorovich problem. Substituting

\bar{ϕ}

into Eq. (3), we can find w is a solution of the elliptic MA equation

ρ_{1} (\nabla w (ξ)) \det \nabla^{2} w (ξ) = ρ_{0} (ξ)

Here, Eq. (5) is a standard MA equation [18].

2.2. A numerical approach of solving the L²Monge-Kantorovich problem

There are two main approaches that can be used to find a numerical solution for Eq. (5). The first one is based on a direct solution of the MA equation [19,20]. The difficulty faced by this approach lies in the discretization of the MA equation. The second approach converts the LMK problem into an evolution problem relying on the gradient of time. For this method, the computational cost may be relevant [15,21]. A numerical method has been presented in [13] for computing the numerical solution of a kind of elliptic MA equation which is a highly nonlinear MA equation in general form [18]. Based on our previous work, an efficient approach will be introduced here for computing the numerical solution of Eq. (5) in order to obtain an appropriate initial design for the MA method.

Instead of solving Eq. (5) directly, we first establish an equivalent problem of Eq. (5), which is given by

\log [ρ_{1} (\nabla w (ξ)) \det \nabla^{2} w (ξ) - ρ_{0} (ξ) + 1] = 0

Notice that ρ₀(ξ) is a normalized density function with max(ρ₀) = 1 on Ω₀. Then, ρ₁(η) is normalized to meet the mass conservation. The value of ρ₀(ξ) could be zero somewhere on Ω₀. That is, ρ₀(ξ) can be a non-negative function, not a strictly positive function. This characteristic is very important for the freeform surface illumination design. We will introduce this characteristic in more detail in Section 3. Equation (6) represents the conservation and redistribution of mass during the mass transport. Obviously, all the inner points of Ω₀ should satisfy Eq. (6). Then, the boundary conditions (BC) are defined to ensure that all the points on the boundary ∂Ω₀ of Ω₀ are mapped to the boundary ∂Ω₁ of Ω₁. Assume ∂Ω₁ can be represented mathematically by the expression f(η) = 0, then BC will be defined as

f (\nabla w (ξ)) = 0 ξ \in \partial Ω_{0}

Since Eq. (7) does not predefine the specific position of each point on ∂Ω₁, consequently the boundary points can automatically move on ∂Ω₁ until the prescribed transport is achieved. Finally, we obtain an equivalent problem of the LMK problem that is given by

{\begin{cases} \log [ρ_{1} (\nabla w (ξ)) \det \nabla^{2} w (ξ) - ρ_{0} (ξ) + 1] = 0 ξ \in {\bar{Ω}}_{0} \\ B C : f (\nabla w (ξ)) = 0 ξ \in \partial Ω_{0} \end{cases}

where

Ω_{0} = {\bar{Ω}}_{0} \cup \partial Ω_{0}

. Then, a numerical method which is similar to the one introduced in [13] is used here to find the solution of Eq. (8). We use the discretization scheme to discretize Eq. (8) and get a set of nonlinear equations. Write the nonlinear equations in the form of

F_{2} (Y) = 0

and linearize this nonlinear problem using the Newton’s method. Further, we can establish an iterative scheme for solving the LMK problem

F_{2} (Y_{k}) + F_{2} ’ (Y_{k}) (Y_{k}_{+ 1} - Y_{k}) = 0

where, Y represents the variables of the nonlinear equations, Y_k denotes the solution obtained from the k-th iteration and F’(Y_k) is the Fréchet derivative of F at Y_k. An initial value is needed to initiate the iteration of Eq. (10). In this new approach, the initial value is defined by

\begin{array}{l} Y_{0} = \frac{1}{2} {[{\bar{X}}_{\min} + \frac{x - x_{\min}}{x_{\max} - x_{\min}} ({\bar{X}}_{\max} - {\bar{X}}_{\min})]}^{2} \frac{x_{\max} - x_{\min}}{{\bar{X}}_{\max} - {\bar{X}}_{\min}} + \frac{1}{2} {[{\bar{Y}}_{\min} + \frac{y - y_{\min}}{y_{\max} - y_{\min}} ({\bar{Y}}_{\max} - {\bar{Y}}_{\min})]}^{2} \\ \times \frac{y_{\max} - y_{\min}}{{\bar{Y}}_{\max} - {\bar{Y}}_{\min}} \end{array}

where, Ω₀₂ = [x_min,x_max] × [y_min,y_max]∈R² and

Ω_{12} = [{\bar{X}}_{\min}, {\bar{X}}_{\max}] \times [{\bar{Y}}_{\min}, {\bar{Y}}_{\max}] \in R^{2}

, as shown in Fig. 2.Ω₀₂ and Ω₁₂ are two rectangular domains which are the minimum rectangles containing Ω₀ and Ω₁, respectively. If Ω₀≠Ω₀₂, the density on Ω₀₁ (Here, Ω₀∪Ω₀₁ = Ω₀₂) is set to be zero. In this case, all the points on Ω₀₁ will converge to ∂Ω₁. Equation (11) shows clearly that the target points defined by the initial mapping will uniformly distribute on Ω₁₂. That is, a uniform density distribution can be produced on Ω₁₂ by this initial mapping if the density distribution defined on Ω₀₂ is uniform. In the next section, we will introduce some characteristics of this new approach, and explore the limitations of the LMK theory in illumination design.

Fig. 2 Define the initial value for the LMK problem.

Download Full Size | PDF

3. Characteristics and limitations

3.1. Characteristics of the new approach

In this subsection, we will introduce some characteristics of the approach presented above, which are relevant to the freeform surface illumination design. To illustrate the rapid convergence of the new approach, we study the following problem

{\begin{cases} ρ_{0} (x, y) = \exp [- 2 (x^{2} + y^{2}) / {0.5}^{2}], ξ = (x, y) \in Ω_{0} = {[- 0.5, 0.5]}^{2} \\ ρ_{1} (t_{x}, t_{y}) = I_{0}, η = (t_{x}, t_{y}) \in Ω_{1} = {[- 2, 2]}^{2} \end{cases}

where ρ₁ represents a uniform density distribution on Ω₁ and the value of I₀ is determined by Eq. (2). Assume Tol is the user defined tolerance, and let Tol = 10⁻¹¹ in all the designs given in this section. The iteration will stop after the i-th iteration when |(∑|F_2,_i|-∑|F_2,_i_-1|)|/N<Tol, where N is the number of the discrete points on Ω₀.

We use 60 × 60 grid points to discretize Ω₀. The stopping criterion is met after seven iterations, as shown in Fig. 3(a).Figure 3(b) gives a gray-scale plot of |F₂| on Ω₀. From this figure, we observe that the |F₂| is of order 10⁻⁴ throughout Ω₀. Here, the numerical error caused by the finite difference scheme is the main reason why |F₂| does not reach zero, as will be explained below. It only takes 11.74 seconds to become saturated by using a computer with a 2.67 GHz Intel(R) Core(TM)2 Quad CPU. Obviously, the iteration converges very fast. Bruneton [9] and Feng [10], respectively, have shown the application of Haker’s method [15] and Sulman’s method [21] in illumination design. These two methods are somewhat similar, because they both rely on the gradient of time. For this kind of method, one should be very careful to choose an appropriate time step dt in order to ensure the convergence. As a comparison between this kind of method and the new approach presented above, we use Sulman’s method to solve the same 2D problem given in Eq. (12). The time step dt is set to be 0.00005. Figure 3(c) shows the convergence of this method. It takes 669.49 seconds to meet the stopping criterion after 28743 iterations. Obviously, the new approach is more efficient than the methods relying on the gradient of time. The example given in Eq. (12) is quite simple. For some other more complex designs, it can be predicted that more time will be needed by the methods relying on the gradient of time. Undoubtedly, the initial design with the LMK theory can benefit from the rapid convergence of the new approach.

Fig. 3 Comparison between the new approach and one kind of method which relies on the gradient of time. (a) Convergence of the new approach, (b) a gray-scale plot of |F₂| on Ω₀ and (c) convergence of the Sulman’s method.

Download Full Size | PDF

Some finite difference schemes are used here to approximate the derivatives included in Eq. (8), and the accuracy of approximation is determined by the spacing h₁ = (x_max-x_min)/m and h₂ = (y_max-y_min)/n in x-axis and y-axis, respectively. Figure 4(a) gives the influence of spacing on the LMK problem given in Eq. (12). Here, we assume m = n. It is easily understood that the minimum mean value of |F₂| will decrease with the decrease of spacing. For example, the ordinate value decreases significantly from 0.00249 to 0.00063 when the spacing decreases from 0.05 to 0.025 (corresponding to the value of m from 20 to 40). However, we also observe that the ordinate value only decreases a little bit when the spacing is already very small. For example, the ordinate value only decreases from 0.00028 to 0.00016 when the spacing decreases from 0.0167 to 0.0125 (corresponding to the value of m from 60 to 80). A smaller spacing means more discrete data points. Thus, it is worthy to make a trade-off between the accuracy of approximation and computational time in illumination design.

Fig. 4 (a) Influence of spacing on the LMK problem. (b) Convergence of the new approach.

Download Full Size | PDF

As stated above, all the points contained in the infinitesimal area element will converge to a point, a straight line or a curve when det(∇ϕ(ξ)) = 0. Here we will discuss this characteristic of the new approach in more detail. Let us consider the following problem

{\begin{cases} ρ_{0} (x, y) = {\begin{cases} \exp [- 2 (x^{2} + y^{2}) / {0.5}^{2}], if \sqrt{x^{2} + y^{2}} \leq {0.5}^{2} \\ 0, elsewise \end{cases}, ξ = (x, y) \in Ω_{0} = {[- 0.5, 0.5]}^{2} \\ ρ_{1} (t_{x}, t_{y}) = I_{0}, η = (t_{x}, t_{y}) \in Ω_{1} \end{cases}

Here, Ω₁ is an ellipse with one-half of the major and minor axes of 4mm and 3mm, respectively. Although the density ρ₀ defined on Ω₀ is not a continuous and strictly position function, the iteration still converges quite fast, as shown in Fig. 4(b). Some discrete points on Ω₀ are predefined in Fig. 5(a), and their corresponding points on Ω₁ produced by the initial design are shown in Fig. 5(b).From this figure, we observe the initial design deviates significantly from the target, and the four boundary corner points of Ω₀₂ are quite far from ∂Ω₁. After the first iteration, the mean value of |F₂| decreases from 0.253 to 0.0246, and the distance between each boundary point of Ω₀₂ and the boundary ∂Ω₁ decreases a lot, especially the four boundary corner points, as shown in Fig. 5(c). Besides, we can also find the inner points of Ω₀₁ have almost already been mapped to ∂Ω₁. After the second iteration, the mean value of |F₂| decreases from 0.0246 to 0.00135, and all the points on Ω₀₁, ∂Ω₀₂ and ∂Ω₀ converge to ∂Ω₁. This design proves the rapid convergence of the new approach again. The characteristic of this approach illustrated in Fig. 5 undoubtedly provides high flexibility for imposing boundary condition. And, this characteristic is especially important for some illumination designs where the incident collimated beam does not have a rectangular cross section or the projected intensity distribution of a compact source is only defined on a circular domain.

Fig. 5 (a) Some discrete points predefined on Ω₀₂, and the distribution of these discrete points obtained from (b) the initial design, (c) the first iteration and (d) the second iteration.

Download Full Size | PDF

3.2. The limitations of the LMK theory in freeform surface illumination design

The unique optimal mapping of the LMK problem is curl-free [15]. Does it mean that one can easily obtain good results just by directly applying the LMK theory to the freeform surface illumination design? We will give the answer in this section, and reveal some limitations of the LMK theory in illumination design. An off-axis design and an on-axis design will be presented here, and the architectures of the two designs are, respectively, plotted in Figs. 6(a) and 6(b). Assume that the incident collimated beam with a uniform intensity distribution should be directed to produce a uniform rectangular pattern on the target plane. The design parameters are given in Table 1.

Fig. 6 Collimated beam shaping with (a) freeform reflector and (b) freeform lens. (c) The beam shaping of a point source. l denotes the distance between the source and the target plane.

Download Full Size | PDF

Table 1. Design parameters (unit: millimeter).

View Table | View all tables in this article

Let us first consider the off-axis design and assume l = 10 mm. At first, we employ the new approach to compute the optimal mapping, as given in Fig. 7(a).According to this mapping, we utilize the PDE method to design the freeform reflector, and obtain the actual mapping produced by the reflector, as shown in Fig. 7(b). Obviously, the actual one deviates significantly from the target one. And, the irradiance distribution produced by the reflector also deviates significantly from the target, as shown in Fig. 8(a).It is clear that the direct application of the LMK theory in this design cannot produce a good result. Although the optimal mapping shown in Fig. 7(a) is curl-free, it may not be curl-free any more in this reflector design due to the optical deflection, and the normal vector field determined by the optimal mapping does not satisfy the integrability condition [8]. As a comparison, we use the MA method to design the reflector, with employing PDE method and the variable separation mapping given in Fig. 7(a) in the initial design of the MA method, and obtain a quite good illumination pattern shown in Fig. 8(b). Since the MA method can automatically satisfy the integrability condition, the mapping obtained from the MA method shown in Fig. 7(c), of course, can be considered to be the optimal mapping of the reflector design. From Figs. 7(a) and 7(c), we can find the two optimal mappings are very different, and the designs respectively obtained from the LMK and the MA method are also quite different.

Fig. 7 (a) The optimal mapping of the LMK problem, and the mappings (b) obtained from the reflector designed by the LMK and (c) the reflector designed by the MA method.

Download Full Size | PDF

Fig. 8 The illumination patterns obtained from (a) the LMK and (b) the MA method. (c) The difference between the intercept point T_ai of the i-th ray and its target point T_ti.

Download Full Size | PDF

We use Dist_i to represent the difference between the intercept point T_ai (tx_i, ty_i, tz_i,) of the i-th incident ray on the target plane and its target point T_ti (Tx_i, Ty_i, Tz_i), which is defined by

D i s t_{i} = \sqrt{{(t x_{i} - T x_{i})}^{2} + {(t y_{i} - T y_{i})}^{2} + {(t z_{i} - T z_{i})}^{2}}

There are usually two normal vector fields that we can obtain when we design the freeform surface with the PDE method or some other methods which rely on a predefined mapping: one is defined by the target mapping and the other one is defined by the smooth freeform surface which is represented by a B-spline surface passing through all the discrete data points. In Fig. 8(c), P_i is an arbitrary point on the freeform surface, T_ti is the target point of point P_i defined by the optimal mapping shown in Fig. 7(a), and N_ti is the target normal vector at point P_i which is determined by P_i and T_ti. N_ai is the actual normal vector at point P_i of the reflector represented by a B-spline surface, and T_ai is the actual intercept point of P_i on the target plane. If the normal vector field defined by the optimal mapping satisfies the integrability condition, the two normal vector fields will be the same regardless of the errors caused by finite difference scheme, data interpolation and construction of the surface. That is, the two normal vectors N_ti and N_ai will equal, T_ti = T_ai and Dist_i = 0. Thus, the value of Dist_i can reflect the integrability of the mapping. From the actual mapping given in Fig. 7(b), we observe the maximum Dist_i is 14.6816 mm, which is approximately 16.4 percent of the length of the diagonal of the rectangular pattern. Obviously, there is much difference between the two normal vector fields. Since the MA method can automatically satisfy the integrability condition, the maximum Dist_i in Fig. 7(c) is very small which is only 0.0425 mm, and only 0.0475 percent of the length of the diagonal.

Change l from 10mm to 100mm, and we can get a relationship between the maximum Dist_i and the lighting distance, as depicted in Fig. 9.Figure 9(a) gives the relationship obtained from the MA method. From this figure, we observe the maximum Dist_i decreases a little bit at the beginning when l increases from 10mm to 20mm. It is easily understood that the profile of the freeform reflector will change more gently with a larger value of l, and the errors caused by finite difference scheme and construction of the surface will be reduced. However, the errors cannot be eliminated even with a larger lighting distance. With further increase of the lighting distance, for example, l increases from 20mm to 100mm, the maximum Dist_i starts to fluctuate irregularly and cannot be further reduced. Figure 9(b) gives the relationship obtained from the LMK. Here, we can find that the distance between the intercept point T_ai and the target point T_ti decreases significantly with the increase of lighting distance. This reflects that the difference between the two normal vector fields is reduced significantly with the increase of lighting distance. Even so, it does not mean the difference can be eliminated. We also change the lighting distance to 1000mm, and find the maximum Dist_i = 0.6766mm, which is much smaller than its initial value of 14.6816mm, but is still larger than the initial value of 0.0425mm obtained from the MA method at l = 10mm. Regardless of those errors, theoretically the difference between the two normal vector fields can be eliminated only and if only the target plane is put at an infinite light distance. The ray-targeting technique which aims to reduce the difference between the actual mapping and the target mapping is usually used to approach a predefined illumination [22]. According to the analyses presented above, one can find the ray-targeting technique may produce some desirable results only when the actual design is close to the target.

Fig. 9 The relationship between the maximum Dist_i and the lighting distance for (a) the MA method and (b) the LMK.

Download Full Size | PDF

Next, let us switch to the on-axis design. By the new approach, the optimal mapping of the LMK problem is obtained, as shown in Fig. 10(a).Assume the refractive index of the freeform lens is 1.59. According to the optimal mapping, we also use the PDE method to design a freeform lens with the lighting distance of 50mm. The mapping and the illumination pattern produced by the freeform lens are, respectively, given in Figs. 10(b) and 11(a).Here, the maximum Dist_i = 4.7716mm. It is clear that the optimal mapping shown in Fig. 10(a) is not curl-free in this lens design due to the refraction of the lens, and there is much difference between the two normal vector fields. The results given in Figs. 10(c) and 11(b) are obtained from the MA method. Obviously, the illumination pattern shown in Fig. 11(b) is quite good. The maximum Dist_i = 0.1148mm which is also much smaller than that obtained from the LMK. Similarly, the mapping obtained from the MA method is also quite different from the optimal mapping of the LMK problem in this on-axis design.

Fig. 10 (a) The optimal mapping of the LMK problem, and the mappings (b) obtained from the lens designed by the LMK and (c) the lens designed by the MA method.

Download Full Size | PDF

Fig. 11 The illumination patterns obtained from (a) the LMK and (b) the MA method. (c) Design flow charts of SMA and LMA methods.

Download Full Size | PDF

Then, we change l from 50mm to 150mm, and obtain a relationship between the maximum Dist_i and the lighting distance, as depicted in Fig. 12.Figure 12(a) gives the relationship obtained from the MA method. Being similar to what we have observed from Fig. 9(a), the maximum Dist_i decreases a little bit at the beginning when l increases from 50mm to 90mm, and starts to fluctuate irregularly with further increase of the lighting distance. Figure 12(b) gives the curve obtained from the LMK. The maximum Dist_i decreases significantly with the increase of lighting distance. This also reflects that the difference between the two normal vector fields is reduced significantly with the increase of lighting distance. Further change the lighting distance to 1000mm, and we get the maximum Dist_i which equals 0.6293mm. Obviously, the maximum Dist_i is much smaller than its initial value of 4.7716mm, however, is still larger than the initial value of 0.1148 mm obtained from the MA method at l = 50mm. There is still a little difference between the two normal vector fields even at l = 1000mm. Similarly, the difference between the two normal vector fields can be eliminated only and if only the target plane is put at an infinite light distance regardless of the errors caused by finite difference scheme, data interpolation and construction of the freeform surface.

Fig. 12 The relationship between the maximum Dist_i and the lighting distance for (a) the MA method and (b) the LMK.

Download Full Size | PDF

According to the detailed analyses of the off-axis and on-axis designs presented above, one can find that the direct application of the LMK theory in illumination design may not be able to produce good results. As introduced in Section 1, the convergence of MA method is strongly determined by the initial design. One common way to obtain the initial design is to utilize a variable separation mapping to construct the freeform surface [5,6,8,13]. In the LMK problem, the optimal mapping is undoubtedly better than a variable separation mapping which is usually not curl-free. Thus, the optimal mapping may be a better option in the initial design of the MA method. In the next section, we will show the significant improvement of the MA method achieved by applying the LMK theory to the initial design.

4. Initial design with the LMK theory

Three design examples will be given here to demonstrate the significant improvement of the MA method achieved by applying the LMK theory to the initial design. In the first example, a freeform reflector is required for directing a collimated beam which has a Gaussian intensity distribution with beam waist of 0.3mm. The target is a uniform elliptical pattern and the architecture of the beam shaping system is given in Fig. 6(a). The other parameters are given in Table 2.For ease of description, we use LMA to represent the MA method which utilizes the optimal mapping of the LMK problem in the initial design, and SMA to represent the MA method which does not utilize the optimal mapping in the initial design. The flow charts of SMA and LMA methods are given in Fig. 11 (c).

Table 2. Design parameters (unit: millimeter).

View Table | View all tables in this article

First, we use the new approach to solve the LMK problem with an initial value defined by Eq. (11), and give the rapid convergence of the LMK problem in Fig. 13(a).After that, we utilize the obtained optimal mapping to design the freeform reflector with PDE method which is used for ray targeting. The illumination pattern produced by the reflector is shown in Fig. 13(b). Since there is still much difference between the two normal vector fields (as analyzed above), the actual illumination pattern deviates from the target. Then, we use this design as the initial value of the MA method, and find the LMA method converges quite fast, as shown in Fig. 13(c). Obviously, the rapid convergence of the LMA method is achieved by using the initial design obtained from the LMK. Besides, a smooth reflector is obtained and the illumination pattern obtained from the LMA is quite good, as given in Fig. 13(d). We also use a variable separation mapping defined by Eq. (11) to find an initial design for the SMA method. The illumination pattern produced by the initial design is given in Fig. 13(e). From this figure, we observe the pattern deviates significantly from that given in Fig. 13(b), and, of course, deviates significantly from the target. With this initial design, we find the iteration cannot converge. More specifically, we encounter complications associated with singularities of the freeform surface after the first iteration and the iteration stops, as shown in Fig. 13(f). From the comparison made between the two initial designs, we can find the LMK provides a better initial design for the LMA method, which is much closer to the final solution. Besides, from Figs. 13(a) and 13(f) we also observe that although the initial mappings used in the LMK problem and the SMA method are the same, the SMA method cannot converge due to the high nonlinearity of the MA equation.

Fig. 13 (a) The convergence of the LMK problem and (b) the illumination pattern obtained from the LMK. (c) The rapid convergence of the LMA method and the obtained smooth reflector, and (d) the illumination pattern obtained from the LMA method. (e) The illumination pattern obtained from a variable separation mapping, (f) the divergence of the SMA method and the singularities of the freeform reflector. Five million rays are traced.

Download Full Size | PDF

In the second example, a freeform lens is required for beam shaping of a point source. The architecture of the beam shaping system is given in Fig. 6(c). The entrance surface of the lens is a spherical surface, and the light source is located at the centre of curvature of the entrance surface. Assume the refractive index of the freeform lens is 1.59, the light source is a Lambertian source, and the maximum collection angle of the lens equals 120°. The first thing we should do to find the optimal mapping of the LMK problem is to get the projected intensity distribution of the point source on a reference plane. Obviously, such a projected intensity distribution is defined on a circular domain. Based on the characteristic of the new approach introduced in Fig. 5, the optimal mapping of the LMK problem is found and the initial design of the freeform lens is obtained. The illumination pattern produced by the initial design is given in Fig. 14(a), and the irradiance distribution along the line x = 0 mm is depicted in Fig. 14(b). These two figures show clearly that the irradiance distribution obtained from the LMK is not quite uniform, and there is much difference between the obtained distribution and the target. With the initial design, we obtain rapid convergence of the LMA method and a quite good approximate to the solution of the freeform lens design, as shown in Figs. 14(c) and 14(d). Then, we use a mapping given in Fig. 14(e) in the initial design. The illumination pattern produced by the initial design is shown in Fig. 14(f). Obviously, the initial design given in Fig. 14(a) is better than that given in Fig. 14(f), and undoubtedly is closer to the target. Moreover, we also find the SMA method cannot iterate at all even at the first iteration due to the total internal reflection taking place on the freeform surface, as shown in Fig. 14(g).

Fig. 14 (a) The illumination pattern obtained from the LMK, and (b) the irradiance curves along the line x = 0mm. (c) The rapid convergence of the LMA method and (d) the illumination pattern obtained from the LMA method. (e) A variable separation mapping which is not the optimal mapping of the LMK problem, (f) the pattern produced by the initial design with five million rays traced and (g) illustration of the total internal reflection taking place on the freeform surface. The unit vector of the incident ray depicted in (g) is (0.8480,0,0.5299) .

Download Full Size | PDF

In the last example, a freeform lens is required for solving the following challenging problem of collimated beam shaping. Figure 6(b) gives the layout of the optical system. The intensity of the incident beam is uniform, and the target irradiance distribution should satisfy

E (x, y) = \frac{C}{3 + \cos [2 π \sqrt{{(x - 0.7)}^{2} + {(y - 0.4)}^{2}}]}

where C is the normalization factor such that Eq. (2) is satisfied. A gray-scale plot of the target distribution is given in Fig. 15(a).Being similar to what we have done in the previous two examples, the first thing is still to use the LMK theory to find an optimal mapping for the initial design. Assume m = n = 80. It takes 60.49 seconds to find the optimal mapping. The illumination pattern produced by the initial design is shown in Fig. 15(b). Since the optimal mapping of the LMK problem is not curl-free in this lens design, the obtained irradiance distribution deviates significantly from the target. Being similar to what Bortz and Shatz have done in [23], we also employ the fractional RMS to quantify the difference between the actual irradiance distribution and the target:

R M S = \sqrt{\frac{1}{M} \sum_{i = 1}^{M} {(\frac{E_{a i} - E_{t i}}{E_{t i}})}^{2}}

where, M is the number of the sample points, E_ti is the target irradiance of the i-th sample point and E_ai is the actual irradiance of the i-th sample point. A smaller value of RMS represents less difference between the actual irradiance distribution and the target. Here, we use 80 × 80 sample points to compute RMS. According to the irradiance distribution, we have RMS = 0.3452 for the initial design. Although there is much difference, we can still observe the illumination pattern produced by the initial design is more or less similar to the target.

Fig. 15 (a) The target illumination pattern, and the illumination patterns obtained from (b) the LMK, (c) the LMA method and (d) the SMA method. Five million rays are traced.

Download Full Size | PDF

Assume that Tol = 10⁻¹¹ and the iteration stops when |(∑|F_1,_i|-∑|F_1,_i_-1|)|/N<Tol. With this initial design obtained from LMK, the LMA method converges quite fast such that the stopping criterion is met after 12 iterations with the computational time of 60.47 seconds, as shown in Fig. 16 (a).The irradiance distribution is shown in Fig. 15(c) and we have RMS = 0.0118. Then, a variable separation mapping defined by Eq. (11) is used here to find the initial design. As mentioned above, the initial design obtained from such a mapping can produce a uniform illumination pattern when the intensity of the incident beam is uniform. With such an initial design, the SMA method also converges and the stopping criterion is met after 18 iterations with the computational time of 94.88 seconds. However, the iteration fluctuates a little bit at the beginning, as shown in Fig. 16(a). The irradiance distribution is shown in Fig. 15(d), and RMS = 0.0117. Since both the LMA method and the SMA method converge to zero in this design regardless of the numerical errors, there is, of course, little difference between the results of the two methods. Here, we need to point out that the SMA method will not converge at l = 8.1mm because we will encounter the complications associated with singularities of the freeform surface, however, the rapid convergence of the LMA method can still be achieved at l = 8.1 mm. Since this case have been studied in the first design example, we do not present the details of this case here. The LMA method converges much faster, and Fig. 16(a) also indicates that the fluctuation can be avoided by using the LMA method.

Fig. 16 Illustration of the convergence of the LMA method and the SMA method at (a) l = 8.2mm and (b) l = 10mm.

Download Full Size | PDF

Change l to 10mm, we find the fluctuation of the SMA is eliminated and both the LMA method and the SMA method converge quite fast, as illustrated in Fig. 16(b). The stopping criterion is met after 12 iterations for the LMA method with the computational time of 63.30 seconds and 13 iterations for the SMA method with the computational time of 73.89 seconds, respectively. Figure 16(a) shows that the stopping criterion is also met after 12 iterations for the LMA method at l = 8.2 mm. This indicates that the convergence of the LMA method is more stable and faster than that of the SMA method.

Although the optimal mapping of the LMK problem may not satisfy the integrability condition in freeform surface illumination design, the optimal mapping is indeed a better option for the initial design of the MA method. The three design examples presented above clearly show that the MA method converges more stably and faster with the application of the LMK theory in the initial design.

5. Conclusion

Due to high nonlinearity of the MA equation, the convergence of the MA method is strongly determined by the initial design. In this paper, we address the initial design of the MA method with the LMK theory. Instead of finding the direct solution of the LMK problem or converting the LMK problem into an evolution problem, we introduce an efficient approach by which the LMK problem is converted into an equivalent problem without relying on the gradient of time. The rapid convergence of this new approach is illustrated and some characteristics of this approach which are relevant to the illumination design are also introduced. Then, some limitations of the LMK theory in illumination design are revealed, and comparisons are made between the LMK theory and the MA method. The optimal mapping of the LMK problem, though, is not curl-free in illumination design, it may still be a good option in the initial design of the MA method. Three design examples, including the beam shaping of collimated beam and point light source, are presented to illustrate the potential benefits of the LMK theory in the initial design of the MA method. Detailed comparisons are made between two different initial designs: one is obtained from the optimal mapping of the LMK problem, and the other one is obtained from the variable separation mapping which is usually used in illumination design. The results clearly show the MA method converges more stably and faster with an initial design obtained from the LMK theory. Moreover, the new approach presented here can also be applied to some other methods of freeform surface illumination design which rely on the initial value.

Acknowledgments

This work was carried out during the tenure of an ERCIM “Alain Bensoussan” Fellowship Programme. The research leading to the results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement n° 246016. Authors thank the European Commission (SMETHODS: FP7-ICT-2009-7 Grant Agreement No. 288526, NGCPV: FP7-ENERGY.2011.1.1 Grant Agreement No. 283798), the Spanish Ministries (ENGINEERING METAMATERIALS: CSD2008-00066, SEM: TSI-020302-2010-65 SUPERRESOLUCION: TEC2011-24019, SIGMAMODULOS: IPT-2011-1441-920000, PMEL: IPT-2011-1212-920000), and UPM (Q090935C59) for the support given to the research activity of the UPM-Optics Engineering Group, making the present work possible.

References and links

1. R. Winston, J. C. Miñano, P. Benítez, N. Shatz, and J. C. Bortz, Nonimaging Optics (Elsevier, 2005).

2. P. Benítez, J. C. Miñano, J. Blen, R. Mohedano, J. Chaves, O. Dross, M. Hernández, and W. Falicoff, “Simultaneous multiple surface optical design method in three dimensions,” Opt. Eng. 43(7), 1489–1502 (2004). [CrossRef]

3. T. L. R. Davenport, T. A. Hough, and W. J. Cassarly, “Optimization for illumination systems: the next level of design,” Proc. SPIE 5456, 81–90 (2004). [CrossRef]

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

5. L. Wang, K. Y. Qian, and Y. Luo, “Discontinuous free-form lens design for prescribed irradiance,” Appl. Opt. 46(18), 3716–3723 (2007). [CrossRef] [PubMed]

6. Y. Ding, X. Liu, Z. R. Zheng, and P. F. Gu, “Freeform LED lens for uniform illumination,” Opt. Express 16(17), 12958–12966 (2008). [CrossRef] [PubMed]

7. V. I. Oliker, “Mathematical aspects of design of beam shaping surfaces in geometrical optics,” in Trends in Nonlinear Analysis, M. Kirkilionis, S. Krömker, R. Rannacher and F. Tomi, eds. (Springer, 2000), pp. 193–224.

8. F. R. Fournier, W. J. Cassarly, and J. P. Rolland, “Fast freeform reflector generation using source-target maps,” Opt. Express 18(5), 5295–5304 (2010). [CrossRef] [PubMed]

9. 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] [PubMed]

10. Z. X. Feng, L. Huang, G. Jin, and M. Gong, “Designing double freeform optical surfaces for controlling both irradiance and wavefront,” Opt. Express 21(23), 28693–28701 (2013). [CrossRef] [PubMed]

11. H. Ries and J. Muschaweck, “Tailored freeform optical surfaces,” J. Opt. Soc. Am. A 19(3), 590–595 (2002). [CrossRef] [PubMed]

12. R. M. Wu, L. Xu, P. Liu, Y. Q. Zhang, Z. R. Zheng, H. F. 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] [PubMed]

13. R. M. Wu, P. Liu, Y. Q. Zhang, Z. R. Zheng, H. F. Li, and X. Liu, “A mathematical model of the single freeform surface design for collimated beam shaping,” Opt. Express 21(18), 20974–20989 (2013). [CrossRef] [PubMed]

14. R. M. Wu, K. Li, P. Liu, Z. R. Zheng, H. F. Li, and X. Liu, “Conceptual design of dedicated road lighting for city park and housing estate,” Appl. Opt. 52(21), 5272–5278 (2013). [CrossRef] [PubMed]

15. S. Haker, L. Zhu, A. Tannenbaum, and S. Angenent, “Optimal mass transport for registration and warping,” Int. J. Comput. Vis. 60(3), 225–240 (2004). [CrossRef]

16. G. L. Delzanno, L. Chacon, J. M. Finn, Y. Chung, and G. Lapenta, “An optimal robust equidistribution method for two-dimensional grid adaptation based on Monge-Kantorovich optimization,” J. Comput. Phys. 227(23), 9841–9864 (2008). [CrossRef]

17. T. Glimm and V. I. Oliker, “Optical design of single reflector systems and the Monge-Kantorovich mass transfer problem,” J. Math. Sci. 117(3), 4096–4108 (2003). [CrossRef]

18. N. S. Trudinger and X.-J. Wang, “The Monge-Ampère equations and its geometric applications,” in Handbook of Geometric Analysis (International Press, 2008), pp. 467–524.

19. E. J. Dean and R. Glowinski, “Numerical methods for fully nonlinear elliptic equations of the Monge Ampere type,” Comput. Methods Appl. Mech. Eng. 195(13–16), 1344–1386 (2006). [CrossRef]

20. B. D. Froese, “A numerical method for the elliptic Monge-Ampere equation with transport boundary conditions,” SIAM J. Sci. Comput. 34(3), A1432–A1459 (2012). [CrossRef]

21. M. M. Sulman, J. F. Williams, and R. D. Russell, “An efficient approach for the numerical solution of Monge-Ampère equation,” Appl. Numer. Math. 61(3), 298–307 (2011). [CrossRef]

22. A. Bruneton, A. Bäuerle, R. Wester, J. Stollenwerk, and P. Loosen, “Limitations of the ray mapping approach in freeform optics design,” Opt. Lett. 38(11), 1945–1947 (2013). [CrossRef] [PubMed]

23. J. Bortz and N. Shatz, “Iterative generalized functional method of nonimaging optical design,” Proc. SPIE 6670, 66700A (2007). [CrossRef]

	x_min	x_max	y_min	y_max	${\bar{X}}_{\min}$	${\bar{X}}_{\max}$	${\bar{Y}}_{\min}$	${\bar{Y}}_{\max}$
Off-axis design	−0.5	0.5	−0.5	0.5	10	50	−40	40
On-axis design	−0.5	0.5	−0.5	0.5	−20	20	−40	40

	x_min	x_max	y_min	y_max	${\bar{X}}_{\min}$	${\bar{X}}_{\max}$	${\bar{Y}}_{\min}$	${\bar{Y}}_{\max}$	l
Example 1	−0.5	0.5	−0.5	0.5	6	8	−3	3	10
Example 2	−0.5	0.5	−0.5	0.5	−5	5	−10	10	10
Example 3	−0.5	0.5	−0.5	0.5	−2	2	−2	2	8.2

	x_min	x_max	y_min	y_max	${\bar{X}}_{\min}$	${\bar{X}}_{\max}$	${\bar{Y}}_{\min}$	${\bar{Y}}_{\max}$
Off-axis design	−0.5	0.5	−0.5	0.5	10	50	−40	40
On-axis design	−0.5	0.5	−0.5	0.5	−20	20	−40	40

	x_min	x_max	y_min	y_max	${\bar{X}}_{\min}$	${\bar{X}}_{\max}$	${\bar{Y}}_{\min}$	${\bar{Y}}_{\max}$	l
Example 1	−0.5	0.5	−0.5	0.5	6	8	−3	3	10
Example 2	−0.5	0.5	−0.5	0.5	−5	5	−10	10	10
Example 3	−0.5	0.5	−0.5	0.5	−2	2	−2	2	8.2

Initial design with L² Monge-Kantorovich theory for the Monge–Ampère equation method in freeform surface illumination design

Abstract

1. Introduction

2. L² Monge-Kantorovich problem and the new approach

2.1. L² Monge-Kantorovich problem

2.2. A numerical approach of solving the L²Monge-Kantorovich problem

3. Characteristics and limitations

3.1. Characteristics of the new approach

3.2. The limitations of the LMK theory in freeform surface illumination design

4. Initial design with the LMK theory

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (16)

Tables (2)

Equations (16)

Optics Express

Abstract

1. Introduction

2. L2 Monge-Kantorovich problem and the new approach

2.1. L2 Monge-Kantorovich problem

2.2. A numerical approach of solving the L2Monge-Kantorovich problem

3. Characteristics and limitations

3.1. Characteristics of the new approach

3.2. The limitations of the LMK theory in freeform surface illumination design

4. Initial design with the LMK theory

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (16)

Tables (2)

Equations (16)

Optics Express

2. L² Monge-Kantorovich problem and the new approach

2.1. L² Monge-Kantorovich problem

2.2. A numerical approach of solving the L²Monge-Kantorovich problem