## Abstract

Computing a source-target map that yields integrable surface normal field is quite challenging for freeform illumination design. Here, we propose a least-squares ray mapping method to calculate a superior ray mapping by iteratively correcting an integrable map to approach the energy conservation and boundary condition. The process is implemented via solving three minimization problems. The first two problems can be figured out pointwise and the third can be converted to two decoupled Poisson equations with Robin boundary conditions. We demonstrate the robustness and high efficiency of the proposed method with several design examples.

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

## Corrections

Shili Wei, Zhengbo Zhu, Zichao Fan, and Donglin Ma, "Least-squares ray mapping method for freeform illumination optics design: erratum," Opt. Express**29**, 11371-11371 (2021)

https://www.osapublishing.org/oe/abstract.cfm?uri=oe-29-7-11371

## 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]:

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).

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:

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}$: 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:

**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

**m**mapping the boundary of $\Omega _{0}$ to the boundary of $\Omega _{1}$. It is necessary to show the following functional:

**t**. Obtaining the matrices $\textbf {P}$ and boundary points

**b**, we can calculate the next predefined map

**m**by minimizing the following functional:

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

#### 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

#### 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 }$:

**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:

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

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**:

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;

#### 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.

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).

#### 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 **A****x**=**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]