Least-squares ray mapping method for freeform illumination optics design

ShiLi Wei; ZhengBo Zhu; ZiChao Fan; DingLin Ma

doi:10.1364/OE.385254

1. Introduction

Designing freeform optics to regulate the spatial energy distribution of a given light source is a fundamental and challenging problem of illumination design. For zero-étendue sources, the design methods can be classified into three main categories [1]: the Monge-Ampère equation method (MA), supporting quadric method (SQM), and the ray mapping method (RM).

Substituting ray-trace equations into the energy conservation equation, the formulation of designing freeform surface to achieve a prescribed irradiance from a point source can be converted to a non-linear partial differential equation of Monge-Ampère type [2–4]. Many methods have been explored to tackle the design problem [4–7], but the theoretical derivation and numerical process are still tedious and complicated. Recently, the iterative wavefront tailoring (IWT) method [8] and the intermediate irradiance transport (IIT) method [9] proposed by Feng can lead to simplified Monge-Ampère equations with simpler mathematical derivation and reduce the design complexity. An alternative approach to calculate freeform surfaces for irradiance control is based on the mathematically rigorous geometric supporting quadric method (SQM) [10,11]. The resulting freeform surface is constructed as a collection of quadric surface patches. For more general cases, the Cartesian oval method [12] is required to solve the illumination design problems.

The most popular method to design freeform illumination optics for zero-étendue sources may be the ray mapping method. Because it is flexible and relatively simple. The ray mapping method divides the illumination design into two steps by firstly computing a source-target ray map and then constructing the freeform optical surfaces to realize the prescribed ray mapping. The most important problem in ray mapping method is to find a ray mapping that the calculated normal field can satisfy the integrability condition [11]:

(1)$$\textbf{N}\cdot(\nabla\times\textbf{N})=0,$$

where N denotes the surface normal vector field. However, it is difficult to find such a ray mapping without solving complex non-linear partial differential equations. The widely used optimal transport mapping with quadric cost function can yield integrable normal field only in paraxial approximation [13–16]. Some optimization methods can be used to improve the illumination performance in non-paraxial cases [13].

Several methods have been developed to address the integrability problem in the ray mapping method. Fournier et al. [11] studied the reflector design by computing an integrable ray mapping with the help of Oliker’s algorithm [10]. Desnijder et al. proposed a method by firstly establishing an optimal transport map [17] and then altering the initial mapping via a symplectic flow to approach the integrability condition [18], which is also demonstrated to be effective in beam shaping system design [19] and freeform lens array design [20]. However, the essence of this method is to eliminate the curl of a rescaled surface normal field N’, which provides a tougher restriction than Eq. (1). This feature may have negative effects on the convergence of the algorithm for different design requirements. Another efficient method is to convert the illumination design to a Monge-Kantorovich mass transfer problem with a special cost function which can be further reduced to a linear assignment problem (LAP) in discrete form [21–24]. Solving the LAP with a special cost function can lead to integrable ray mapping for a special optical system. However, the discretized search space may potentially increase the computational expense for complex irradiance problems, such as picture-generated illumination system design.

In this paper, we propose a least-squares ray mapping (LSRM) method to compute a superior source-target map that the calculated normal field satisfies the integrability condition. This method traces the rays redirected by the freeform surface to calculate an integrable target map and iteratively corrects this map to satisfy energy conservation and boundary condition. We elaborate on the basic idea and design process in Sections 2 and 3. In Section 4, several challeging design examples are achieved using the proposed LSRM method. In Section 5, a brief conclusion is presented.

2. Formulation of the problem and idea

2.1 Formulation of the problem

Zero-étendue beam is defined as a normal system of rays, i.e., there is a family of surfaces normal to the rays (the wavefronts) [25]. Therefore, each ray of a zero-étendue beam can be determined by two parameters. Assume that the source plane is defined in Cartesian coordinates $(u, v)$ at $z=z_{0}$ and let $f$ be the irradiance in a bounded region $\Omega _{0}$. Similarly, we can define the target plane in coordinates $(x, y)$ at $z=z_{1}$ and $g$ is the prescribed irradiance distribution in region $\Omega _{1}$. The input ray bundle can be defined using coordinates $(u, v)$. Let $\textbf {W}(u, v)$ be the input wavefront and $\textbf {I}(u, v)$ be the normalized incident ray vector. The inverse problem requires us to design a freeform surface $\textbf {L}=(L_{x},L_{y},L_{z})$ to redirect the rays $\textbf {I}(u, v)$ to the target plane so that we can obtain a predefined irradiance $g(x, y)$. The sketch is presented in Fig. 1(a).

Fig. 1. (a) Freeform illumination design for zero-étendue sources; (b) the predefined optimal transport target map and the real target map; (c) the simulated irradiance of the designed system using optimal transport mapping.

Download Full Size | PDF

In ray mapping method, we need to find a map $\textbf {m}:\Omega _{0}\rightarrow \Omega _{1}$ redistributing the energy density $f$ to $g$. Using the energy conservation for one infinitesimal light tube, we can obtain:

(2)$$g(\textbf{m}(u, v))\textrm{det}(D\textbf{m}(u, v))=f(u, v),$$

where $D\textbf {m}$ denotes the Jacobi matrix of m. Furthermore, the boundary rays of $\Omega _{0}$ should be mapped to the boundary of $\Omega _{1}$:

(3)$$\textbf{m}(\partial\Omega_{0})=\partial\Omega_{1}.$$

The most critical and intractable problem is to find a mapping m that can simutaneously satisfy Eqs. (1)–(3).

2.2 How to compute an integrable source-target map

We first consider a design example redirecting rays from a Lambertian point source to form a uniform distribution with non-convex boundary. The point source is located at the origin point. The source produces an irradiance distribution proportional to $z_{0}^{4}/(u^{2}+v^{2}+z_{0}^{2})^{2}$ on the source plane $\Omega _{0}$ at $z=z_{0}=1$ mm. The source domain is defined as: $\Omega _{0}=\{(u, v) \big | |u| \leq 1 \textrm {mm}, |v| \leq 1 \textrm {mm}\}$. The target boundary is defined by $\rho (\theta )=\rho _{\textrm {max}}(1+0.2\textrm {cos}(3\theta )), 0\leq \theta \leq 2\pi ,$ where $\rho$ is the distance to the origin at plane $z=z_{1}$ and where $\theta$ is the counterclockwise angle with respect to the $x$-axis. We set $\rho _{\textrm {max}}$=80 mm and $z_{1}$=60 mm. To achieve the design requirements, a common method is to employ the optimal transport mapping. We set a 151$\times$151 regular source grid and calculate an optimal transport target map shown as the red grid in Fig. 1(b). We employ a least-square surface construction method proposed by Feng [16] to compute the freeform surface $\textbf {L}$. Other surface construction methods can also be used in this paper. However, the optimal transport map does not yield an integrable surface normal field in the non-paraxial case which means that the calculated surface can not redirect the rays to the correct position on the target plane. If we obtain the surface points, we can calculate the surface normal field:

(4)$$\textbf{N}=\frac{\frac{\partial\textbf{L}}{\partial u}\times\frac{\partial\textbf{L}}{\partial v}}{|\frac{\partial\textbf{L}}{\partial u}\times\frac{\partial\textbf{L}}{\partial v}|}.$$

The refracted rays can be simply computed via Snell’s law. Therefore, we can trace the rays from the source to calculate a real target map t shown as the blue grid in Fig. 1(b). As presented in Fig. 1(b), the blue grid t does not totally overlap the red grid m. It seems that the real map does not deviate from the predefined map too much, but the irradiance distribution is seriously degraded from the desired target distribution as shown in Fig. 1(c).

Through the above process, we can obtain two maps. The predefined map m$(m_1, m_2)$ satisfies the energy conservation Eq. (2) and boundary condition Eq. (3) but does not meet the integrability condition Eq. (1). The real target map t$(X,Y)$ satisfies integrability condition Eq. (1) but does not meet Eqs. (2) and (3). Two ways are promising to achieve a favorable map satisfying Eqs. (1)–(3). One way is to alter the predefined map m to approach the integrablity condition without changing the satisfaction of Eqs. (2) and 3). And this process can be achieved via an area-preserving transformation which is presented by Desnijder [18].

In this paper, we propose another approach to compute a superior ray mapping satisfying Eqs. (1)–(3) by focusing on the real map t. The problem is whether we can change the integrable real map t to satisfy the energy conservation Eq. (2) and boundary condition Eq. (3). We show in the paper that this requirement can be achieved using a least-squares algorithm [26] with some essential modifications. Note that we do not need to calculate an optimal transport mapping as an initial mapping in this method. We can just simply establish a regular grid and the algorithm can also converge well as presented in Sections 3 and 4.

3. Least-squares ray mapping method

The original least-squares algorithm is presented in [26] which is used to solve the optimal transport problem. In this paper, we modify the algorithm and combine it into our design framework to calculate a superior ray mapping. The method starts with an initial mapping $\textbf {m}^{0}$. The algorithm is not sensitive to the starting point, which means that the initial mapping does not need to be defined by some particular calculation. We can just establish a rectangular grid that has a similar size to the illumination region as presented in Section 4.1.

Assume that $\textbf {m}=\textbf {m}^n$ is given, and we can construct the freeform surface and compute a real map $\textbf {t}^n$. The goal is to correct the map $\textbf {t}^n$ to satisfy Eqs. (2) and (3). Therefore, we require the Jacobi matrix of $\textbf {t}$ to be given by

(5)$$D\textbf{t}= \begin{pmatrix} \frac{\partial{X}}{\partial{u}} & \frac{\partial{X}}{\partial{v}}\\ \frac{\partial{Y}}{\partial{u}} & \frac{\partial{Y}}{\partial{v}} \end{pmatrix},$$

equals a positive semidefinite matrix $\textbf {P}(u, v)$ satisfying:

(6)$$\textrm{det}(\textbf{P}(u, v))=\frac{f(u, v)}{g(\textbf{m}(u, v))},$$

for all $(u, v)\in \Omega _{0}$. Note that we discard the requirement of the matrix $\textbf {P}$ to be symmetric in [26] because the map can not be regarded as a gradient of a convex function any more. We need to find matrices $\textbf {P}$ that satisfies Eq. (6) and minimizes the following functional:

(7)$$J_{1}(\textbf{t},\textbf{P})=\frac{1}{2}\iint_{\Omega_{0}}\left\|D\textbf{t}-\textbf{P}\right\|^2dudv,$$

where $\left \|A\right \|^2$ is defined as $\sum a_{i, j}^2$. Besides, we require m mapping the boundary of $\Omega _{0}$ to the boundary of $\Omega _{1}$. It is necessary to show the following functional:

(8)$$J_{B}(\textbf{t},\textbf{b})=\frac{1}{2}\oint_{\partial\Omega_{0}}\left|D\textbf{t}-\textbf{b}\right|^2ds.$$

The functional Eq. (8) is minimized to find the boundary points which are closest to the current real map t. Obtaining the matrices $\textbf {P}$ and boundary points b, we can calculate the next predefined map m by minimizing the following functional:

(9)$$J(\textbf{m},\textbf{P},\textbf{b})=(1-\alpha)J_{B}(\textbf{m},\textbf{b})+\alpha J_{1}(\textbf{m},\textbf{P}),$$

where the parameter $\alpha$ controls the weight of $J_1$ to $J_B$.

We start with an initial mapping $\textbf {m}^{0}$ and perform the iteration subsequently:

(10)$$\begin{aligned} & \textbf{b}^{n+1}=\arg\min_{\textbf{b}}J_{B}(\textbf{t}^{n},\textbf{b}),\\ & \textbf{P}^{n+1}=\arg\min_{\textbf{P}}J_{1}(\textbf{t}^{n},\textbf{P}),\\ & \textbf{m}^{n+1}=\arg\min_{\textbf{m}}J(\textbf{m},\textbf{P},\textbf{b}). \end{aligned}$$

When we obtain the updated predefined map $\textbf {m}^{n+1}$, we can use it to calculate the approximate freeform surface and then compute the real target map $\textbf {t}^{n+1}$. The iteration will loop until $\textbf {m}$ getting close to $\textbf {t}$. The design process is shown in Fig. 2.

Fig. 2. The iteration flow chart.

Download Full Size | PDF

3.1 Minimization procedure for boundary points $\textbf {b}$

Assume that $\textbf {t}=\textbf {t}^n$ is given, and we minimize the functional $J_{B}$ to find the closest boundary points. This process can be performed pointwise. We discretize the boundary of $\Omega _{1}$ by defining many points on $\partial \Omega _{1}$. We connect adjacent points by line segments and determine the closest point to $\textbf {t}_{i, j}\in \partial \Omega _{1}$ on all line segments. And we can calculate the nearest point $\textbf {b}_{i, j}$ to all points $\textbf {t}_{i, j}\in \partial \Omega _{1}$. This process is the same as Section 3.1 in [26].

3.2 Minimization procedure for matrices $\textbf {P}$

In this section, we perform the minimization procedure for matrices $\textbf {P}$. The mathematical treatment is analogic to Section 3.2 in [26]. The main difference is that we do not require the matrices $\textbf {P}$ to be symmetric. We introduce the notation

(11)$$D\textbf{t}= \begin{pmatrix} d_{11} & d_{12}\\ d_{21} & d_{22} \end{pmatrix}, \quad\textbf{P}= \begin{pmatrix} p_{11} & p_{12}\\ p_{21} & p_{22} \end{pmatrix},$$

with $d_{11}=\delta _{u}X$, $d_{12}=\delta _{v}X$, $d_{21}=\delta _{u}Y$, $d_{22}=\delta _{v}Y$, where $\delta _{u}$ and $\delta _{v}$ are the central difference approximations of $\partial /\partial u$ and $\partial /\partial v$. Next, we introduce the function:

(12)$$H(p_{11}, p_{12}, p_{21}, p_{22})=\left\|D\textbf{t}-\textbf{P}\right\|^2=\sum_{k, l}(d_{kl}-p_{kl})^2.$$

For every point $(u_{i, j}, v_{i, j})$ in $\Omega _{0}$, the matrix $\textbf {P}_{i, j}$ should satisfy Eq. (6) and minimize Eq. (12). The possible minimizers are critical points of the Lagrangian function $L=L(p_{11}, p_{12}, p_{21}, p_{22};\lambda )$ defined by

(13)$$\begin{aligned} L(p_{11}, p_{12}, p_{21}, p_{22};\lambda) & =H(p_{11}, p_{12}, p_{21}, p_{22})+\lambda(\textrm{det}(\textbf{P})-\frac{f}{g})\\ & =(d_{11}-p_{11})^2+(d_{12}-p_{12})^2+(d_{21}-p_{21})^2+(d_{22}-p_{22})^2+\lambda(\textrm{det}(\textbf{P})-\frac{f}{g}), \end{aligned}$$

where $\lambda$ is the Lagrange multiplier. Setting all partial derivatives of $L$ to be 0, we obtain the following algebraic system:

(14a)$$p_{11}+\lambda p_{22}=d_{11},$$

(14b)$$p_{12}-\lambda p_{21}=d_{12},$$

(14c)$$p_{21}-\lambda p_{12}=d_{21},$$

(14d)$$p_{22}+\lambda p_{11}=d_{22},$$

(14e)$$p_{11}p_{22}-p_{12}p_{21}=\frac{f}{g},$$

From Eq. (14) we can obtain the value of each element in matrix $\textbf {P}$:

(15)$$p_{11}=\frac{\lambda d_{22}-d_{11}}{\lambda^2-1}, p_{12}=\frac{d_{12}+\lambda d_{21}}{1-\lambda^2}, p_{21}=\frac{\lambda d_{12}+d_{21}}{1-\lambda^2}, p_{22}=\frac{\lambda d_{11}-d_{22}}{\lambda^2-1}.$$

Note that we do not consider the case that $\lambda =\pm 1$ in this paper. $\lambda =-1$ occurs when $d_{11}+d_{22}=0$, which will lead impractical crossing grid line. The case that $\lambda =1$ ($d_{11}=d_{22}, d_{12}=-d_{21}$) almost never happens in our design framework. The case that $d_{11}=d_{22}, d_{12}=-d_{21}$ can be treated with the same process as [26]. And the algorithm can quickly revise the mapping to the correct position in the next iterative steps. Substituting Eq. (15) to Eq. (14e) leads to a quartic equation:

(16)$$a_{4}\lambda^4+a_{2}\lambda^2+a_{1}\lambda+a_{0}=0.$$

Solving the quartic equation, we can obtain four possible minimizers for Eq. (12). Substituting the solutions to Eq. (12) and selecting the correct minimizer, the matrix $\textbf {P}_{i, j}$ can be computed through Eq. (15). Performing the above process pointwise, we can obtian the matrices $\textbf {P}$.

3.3 Minimization procedure for target map m

In this section, we calculate the mapping m to minimize the functional Eq. (9). This step cannot be implemented pointwise. We calculate the first variation $\delta J(\textbf {m},\textbf {P},\textbf {b})(\eta )$ with respect to m for $\boldsymbol {\eta }$:

(17)$$\begin{aligned} \delta J(\textbf{m},\textbf{P},\textbf{b})(\eta) & =\lim_{\epsilon\to0}\frac{1}{\epsilon}(J(\textbf{m}+\epsilon\boldsymbol{\eta}, \textbf{P},\textbf{b})- J(\textbf{m},\textbf{P},\textbf{b}))\\ & =\lim_{\epsilon\to0}[\frac{\alpha}{2}\iint_{\Omega_{0}}2(D\textbf{m}-\textbf{P}):D\boldsymbol{\eta}+\epsilon\left\|D\boldsymbol{\eta}\right\|^2dudv\\ & +\frac{1-\alpha}{2}\oint_{\partial\Omega_{0}}2(\textbf{m}-\textbf{b})\cdot\boldsymbol{\eta}+\epsilon|\boldsymbol{\eta}|^2ds]\\ & =\alpha\iint_{\Omega_{0}}(D\textbf{m}-\textbf{P}):D\boldsymbol{\eta}dudv+(1-\alpha)\oint_{\partial\Omega_{0}}(\textbf{m}-\textbf{b})\cdot\boldsymbol{\eta}ds. \end{aligned}$$

The notation "A:B" denotes the Fröbenius inner product of the matrices $\textbf {A}=(a_{i, j})$ and $\textbf {B}=(b_{i, j})$, defined by $\textbf {A}:\textbf {B}=\sum {a_{i, j}b_{i, j}}$. The minimizer of Eq. (17) is given by $\delta J(\textbf {m},\textbf {P},\textbf {b})(\eta )=0$. Applying Gaussian’s theorem and fundamental lemma of the calculus of variations [27], we obtain two decoupled Poisson equations with Robin boundary conditions:

(18a)$$\Delta m_1=\nabla\cdot\textbf{p}_1,\qquad(u, v)\in\Omega_0,$$

(18b)$$(1-\alpha)m_1+\alpha\nabla m_1\cdot\hat{\textbf{n}}=(1-\alpha)b_x+\alpha\textbf{p}_1\cdot\hat{\textbf{n}},\quad(u, v)\in\partial\Omega_0,$$

(18c)$$\Delta m_2=\nabla\cdot\textbf{p}_2,\qquad(u, v)\in\Omega_0,$$

(18d)$$(1-\alpha)m_2+\alpha\nabla m_2\cdot\hat{\textbf{n}}=(1-\alpha)b_y+\alpha\textbf{p}_2\cdot\hat{\textbf{n}},\quad(u, v)\in\partial\Omega_0,$$

where $\textbf {p}_1=(p_{11}, p_{12})^T$, $\textbf {p}_2=(p_{21}, p_{22})^T$ and $\hat {\textbf {n}}$ is the outward unit normal vector on $\partial \Omega _{0}$. $b_x$ and $b_y$ denote the Cartesian coordinates of boundary points calculated in Section 3.1. We solve the Poisson boundary value problems by standard second-order differences for the first and second-order derivatives. The derivatives of $p_{11}$, $p_{12}$, $p_{21}$, and $p_{22}$ are approximated using the same finite difference schemes. When calculating the derivatives of points on the boundary, we use the one-sided second-order schemes. The Poisson equation can be transformed to a system of linear equations. Solving the group of linear equations, we can obtain the coordinates of updated mapping m.

Performing the iteration Eq. (10) with the procedures demonstrated in Sections 3.1 to 3.3, we can alter an initial integrable target map to approach the energy conservation Eq. (2) and the boundary condition Eq. (3), which will yield a superior ray mapping satisfying Eqs. (1)–(3) for freeform illumination design. We will present several challenging design examples in Section 4 to demonstrate the robustness and efficiency of the proposed LSRM method.

4. Design examples

4.1 Iteration process of LSRM

We first focus on the design examples demonstrated in Section 2.2. The source produces an irradiance distribution proportional to $z_{0}^{4}/(u^{2}+v^{2}+z_{0}^{2})^{2}$ on the source plane $\Omega _{0}$ at $z=z_{0}=1$ mm. The source domain is defined as: $\Omega _{0}=\{(u, v) \big | |u| \leq 1 \textrm {mm}, |v| \leq 1 \textrm {mm}\}$. The target boundary is defined by $\rho (\theta )=\rho _{\textrm {max}}(1+0.2\textrm {cos}(3\theta )), 0\leq \theta \leq 2\pi ,$ where $\rho$ is the distance to the origin at plane $z=z_{1}$ and $\theta$ is the counterclockwise angle with respect to the $x$-axis. We set $\rho _{\textrm {max}}$=80 mm and $z_{1}$=60 mm. The lens material is polymethyl methacrylate (PMMA). We set a 151$\times$151 uniform grid on the source plane and $\alpha$=0.05. There is no need to calculate an optimal transport map as an initial mapping. We can just start with a rectangular grid with similar size of illumination region and the algorithm converges well as shown in Fig. 3(a). Other reasonable initial mapping can also be used to implement the algorithm. The choice of starting points just influences the convergence time of the algorithm but do not change the final result.

Fig. 3. (a) The convergence process of predefined map m and real map t; (b) the RMS deviation between predefined map m and real map t as well as computation time; (c) simulated target irradiance distribution produced by the resulted freeform surface.

Download Full Size | PDF

We perform the iteration 40 times and the real map t is getting close to the predefined map m as presented in Fig. 3(a). We use a root-mean-square (RMS) value to evaluate the deviation of map m and t:

(19)$$\delta_{1}=\sqrt{\frac{1}{mn}(\sum_{i}\sum_{j}|\textbf{m}_{i, j}-\textbf{t}_{i, j}|^2)},$$

where $m$ and $n$ denote the maximum values of the grid number in both directions. As presented in Fig. 3(b), the deviation is reduced from 1.62 mm in the first iteration to 0.015 mm after 40 iteration processes. The computation time of 40 iterations is 123s. And the deviation value can be further reduced to $7.4\times 10^{-3}$ mm after 100 iterations. We can further obtain a smaller deviation value via setting a finer grid or performing more iterations. The simulated irradiance distribution is shown in Fig. 3(c), which presents a significant improvement compared to the optimal transport map shown in Fig. 1(c).

The parameter $\alpha$ controls the weight of $J_{1}$ to $J_{B}$. It influences the algorithm’s rate of convergence to approach Eqs. (1)–(3). The integrability condition can be evaluated by the root-mean-square deviation (RMSD) between predefined map m and real map t, which is defined by Eq. (19). The energy conservation Eq. (2) can be evaluated by an RSMD value defined by;

(20)$$\delta_{2}=\sqrt{\frac{1}{mn}[\sum_{i}\sum_{j}(\textrm{det}(D\textbf{t}_{i, j})g(\textbf{t}_{i, j})-f(u_{i, j}, v_{i, j}))^2]},$$

Similary, the boundary condition Eq. (3) can be calculated by an RSMD value defined by;

(21)$$\delta_{3}=\sqrt{\frac{1}{2m+2n-4}[\sum_{i=1, m}\sum_{j=1, n}(\textbf{t}_{i, j}-\textbf{b}_{i, j})^2]},$$

We perform the iteration with different values of $\alpha$ and test the convergence rate of the proposed algorithm as presented in Fig. 4. As shown in Figs. 4(b) and 4(c), if we set a small value of $\alpha$, the algorithm can converge quickly to approach the boundary condition but do not perform well in the convergence speed of energy conservation. Conversely, the algorithm approaches the energy conservation fast if we set a relatively large $\alpha$ value. But the convergence of boundary condition Eq. (3) will cost more iteration processes. Based on our experience, selecting the value of $\alpha$ to be between 0.03 and 0.2 will be appropriate.

Fig. 4. (a) The RMSD from predefined map m to real map t; (b) the RMSD of energy conservation; (c) the RMSD of boundary points.

Download Full Size | PDF

4.2 More design examples

The next two design examples are aimed at producing a uniform rectangular illumination and generating a letter "S" on a non-convex boundary region. The rectangular region is defined as $\Omega _{1}=\{(x, y) \big | |x| \leq 80 \textrm {mm}, |y| \leq 40 \textrm {mm}\}$ at $z_1=100$ mm. The non-convex boundary is defined by $\rho (\theta )=\rho _{\textrm {max}}(1+0.15\textrm {cos}(6\theta ))$, where $\rho _{\textrm {max}}=80 \textrm {mm}$ and $z_1$ is set to be 60 mm. The irradiance ratio is set to be 3:1 and $\alpha$ is set to be 0.05. The source domain and irradiance $f$ are the same as the first design example. We use a 201$\times$201 grid and implement the LSRM method demonstrated in Section 3. As shown in Figs. 5(a) and 5(c), the designed illumination system is beyond the paraxial approximations. The resulted source-target maps are presented in Figs. 5(b) and 5(d). The simulated irradiance distributions are shown in Figs. 5(e) and 5(f). The computation time of one iteration is 0.75s for a 101$\times$101 grid and 21s for a 201$\times$201 grid. We can first use a 101$\times$101 grid to perform the iterative step dozens of times, then interpolate a 201$\times$201 grid and perform several further iterations. The whole design process is completed in several minutes.

Fig. 5. (a,c) The designed illumination optics system; (b,d) the result source-target maps; (e,f) the simulated irradiance distributions for uniform rectangular pattern and letter pattern in a non-convex illuminated region.

Download Full Size | PDF

To verify the robustness of the proposed method, we design a picture-generated illumination system. The target plane is located at $z_1=100$ mm. The illumination region is a rectangular area defined by $\Omega _{1}=\{(x, y) \big | |x| \leq 62.5 \textrm {mm}, |y| \leq 110 \textrm {mm}\}$. We define an elliptic paraboloid entrance surface with analytical formula $z=z_0(1-x^2/x_0^2-y^2/y_0^2)$. And we set $z_0=15$ mm, $x_0=25$ mm and $y_0=15$ mm. The original picture is presented in Fig. 6(c). The result illumination system and the target map are shown in Figs. 6(a) and 6(b). The simulated irradiance is presented in Fig. 6(d). Notice that we only show a 101$\times$101 grid in Fig. 6(b) for better visualization. The algorithm is implemented in a 401$\times$401 grid. And we can obtain a better result than Fig. 6(d) via defining a finer grid. We present the resulting freeform lens and the Gaussian curvature of the freeform surface in Figs. 6(e) and 6(f). We first start the algorithm with a 101$\times$101 grid and perform the iteration 29 times in less than 25s. Then we interpolate a 201$\times$201 grid and perform iterative steps 5 times. This process cost approximate 133s. In the last, a 401$\times$401 grid is interpolated and the iteration is implemented 3 times which spends 1836s to obtain the final result as presented in Fig. 6(g).

Fig. 6. (a) The designed illumination system; (b) the result source-target map; (c) the original picture; (d) the simulated irradiance distribution; (e) the freeform lens; (f) Gaussian curvature of the freeform surface; (g) the iteration process and computation time.

Download Full Size | PDF

4.3 Discussions

The computation processes of the above design examples are performed on a standard computer with an Intel Core i7 8700k and 32GB RAM. The most time-consuming processes are the surface construction procedure and solving two Poisson equations. Since they all require to solve a large system of linear equations Ax=b. We use a sparse matrix to store the coefficient matrix A and perform function $mldivide$ in Matlab to calculate the surface points and the coordinates of the updated mapping m. Therefore, the time cost increases steeply when adding the number of grid points. We can further use some efficient algorithms to solve the large group of linear equations to reduce the time complexity of this method.

The method considers the geometry of the freeform surface in the calculation of ray mapping by constructing the freeform surface in every iterative step. Therefore, the method can both tackle the illumination design problems in cases that the size of the freeform lens can be neglected or not (i.e. far-field approximation or not).

The current least-squares ray mapping method (LSRM) proposed in this paper is implemented in Cartesian coordinates $(u, v)$, which presents high efficiency in computing a superior ray mapping to satisfy the integrability condition. The method can be further extended to other orthogonal parametric spaces such as polar coordinates $(\rho ,\theta )$ to adapt the light source domain. The possible mathematical treatment can refer to [7,28]. If the source boundary is arbitrary, some other mathematical treatment is required to accomplish the design. And we will consider this point in our future work.

5. Conclusion

We proposed an efficient LSRM method to address the integrability problem in the ray mapping method for freeform illumination design. The method is to iteratively correct an integrable map to approach energy conservation and boundary condition using a modified least-squares algorithm. Benefit from the heuristic processes, the method can calculate a superior ray mapping without solving complex non-linear partial differential equations of Monge-Ampère type. The algorithm is implemented by two pointwise steps and solving two Poisson equations which can be converted to two systems of linear equations. The computation time is also within tolerance as shown in the design examples. We demonstrated the efficiency and robustness of the proposed method by designing freeform lenses to generate complex illuminated patterns.

Funding

National Natural Science Foundation of China (61805088); Fundamental Research Funds for the Central Universities (2019kfyRCPY083, 2019kfyXKJC040).

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. J. Schruben, “Formulation of a reflector-design problem for a lighting fixture,” J. Opt. Soc. Am. 62(12), 1498–1501 (1972). [CrossRef]

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

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

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. C. Bösel and H. Gross, “Single freeform surface design for prescribed input wavefront and target irradiance,” J. Opt. Soc. Am. A 34(9), 1490–1499 (2017). [CrossRef]

7. L. B. Romijn, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Freeform lens design for a point source and far-field target,” J. Opt. Soc. Am. A 36(11), 1926–1939 (2019). [CrossRef]

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

9. Z. Feng, D. Cheng, and Y. Wang, “Transferring freeform lens design into phase retrieval through intermediate irradiance transport,” Opt. Lett. 44(22), 5501–5504 (2019). [CrossRef]

10. 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-Verlag, 2003), pp. 193–224.

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

12. D. Michaelis, P. Schreiber, and A. Bräuer, “Cartesian oval representation of freeform optics in illumination systems,” Opt. Lett. 36(6), 918–920 (2011). [CrossRef]

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

14. C. Bösel and H. Gross, “Ray mapping approach for the efficient design of continuous freeform surfaces,” Opt. Express 24(13), 14271–14282 (2016). [CrossRef]

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

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

17. K. Desnijder, P. Hanselaer, and Y. Meuret, “Flexible design method for freeform lenses with an arbitrary lens contour,” Opt. Lett. 42(24), 5238–5241 (2017). [CrossRef]

18. K. Desnijder, P. Hanselaer, and Y. Meuret, “Ray mapping method for off-axis and non-paraxial freeform illumination lens design,” Opt. Lett. 44(4), 771–774 (2019). [CrossRef]

19. S. Wei, Z. Zhu, Z. Fan, Y. Yan, and D. Ma, “Double freeform surfaces design for beam shaping with non-planar wavefront using an integrable ray mapping method,” Opt. Express 27(19), 26757–26771 (2019). [CrossRef]

20. K. Desnijder, W. Ryckaert, P. Hanselaer, and Y. Meuret, “Luminance spreading freeform lens arrays with accurate intensity control,” Opt. Express 27(23), 32994–33004 (2019). [CrossRef]

21. L. L. Doskolovich, A. A. Mingazov, D. A. Bykov, E. S. Andreev, and E. A. Bezus, “Variational approach to calculation of light field eikonal function for illuminating a prescribed region,” Opt. Express 25(22), 26378–26392 (2017). [CrossRef]

22. 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]

23. D. A. Bykov, L. L. Doskolovich, A. A. Mingazov, E. A. Bezus, and N. L. Kazanskiy, “Linear assignment problem in the design of freeform refractive optical elements generating prescribed irradiance distributions,” Opt. Express 26(21), 27812–27825 (2018). [CrossRef]

24. 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]

25. J. C. Minano, “Application of the conservation of etendue theorem for 2-d subdomains of the phase space in nonimaging concentrators,” Appl. Opt. 23(12), 2021–2025 (1984). [CrossRef]

26. C. Prins, R. Beltman, J. ten Thije Boonkkamp, W. IJzerman, and T. W. Tukker, “A least-squares method for optimal transport using the monge–ampère equation,” SIAM J. Sci. Comput. 37(6), B937–B961 (2015). [CrossRef]

27. R. Courant and D. Hilbert, Methods of Mathematical Physics: Partial Differential Equations (John Wiley & Sons, 2008).

28. R. Beltman, J. ten Thije Boonkkamp, and W. IJzerman, “A least-squares method for the inverse reflector problem in arbitrary orthogonal coordinates,” J. Comput. Phys. 367, 347–373 (2018). [CrossRef]

Least-squares ray mapping method for freeform illumination optics design

Abstract

Corrections

1. Introduction

2. Formulation of the problem and idea

2.1 Formulation of the problem

2.2 How to compute an integrable source-target map

3. Least-squares ray mapping method

3.1 Minimization procedure for boundary points $\textbf {b}$

3.2 Minimization procedure for matrices $\textbf {P}$

3.3 Minimization procedure for target map m

4. Design examples

4.1 Iteration process of LSRM

4.2 More design examples

4.3 Discussions

5. Conclusion

Funding

Disclosures

References

Cited By

Figures (6)

Equations (28)

Optics Express