Open Access
Add to BibSonomy
Abstract
The ray-mapping method has been widely used for designing freeform illumination lenses. However, in non-paraxial or off-axis situations, it remains challenging to obtain an integrable ray-mapping, often requiring a complex iterative correction process for the initial mapping. To address this challenge, we propose an extended ray-mapping method that incorporates differentiable ray-tracing into the design pipeline of the ray-mapping method. This enables accurate surface construction according to ray-mapping and efficient shape correction based on irradiance distribution. The proposed method involves two optimization stages. In the first stage, the freeform surface is preliminarily optimized to closely match the optimal transport mapping. The obtained freeform surface is then further optimized in the second stage to minimize the divergence between the target and simulated irradiance distributions. Additionally, the mean curvature of the freeform surface is also constrained in the second stage to facilitate the fabrication of the final freeform surface. Non-paraxial illumination lenses and off-axis illumination lenses have been designed using the proposed method within ten minutes, and simulations demonstrate that the approach is effective and robust.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Designing freeform optical surfaces that can generate a specified irradiance distribution on the target plane from the given intensity distribution of a light source is one of the fundamental yet challenging problems in non-imaging optics. Three crucial constraints must be considered simultaneously in the freeform lenses design. Firstly, smoothness constraint must be taken into account to ensure manufacturability, which describes the integrable condition of the normal vector $\mathbf {N}$ of the surface [1]:
Secondly, light propagation through freeform surfaces must satisfy Snell’s law. Thirdly, energy conservation between the source and the target planes should be guaranteed, which is expressed as a partial differential equation [2]: where $\mathbf{\phi }: (x_s,y_s)\rightarrow (x_t,y_t)$ is a map between the source and the target planes, while $\nabla \mathbf{\phi }$ represent the Jacobian matrix of $\mathbf{\phi }$. $I_s(x_s, y_s)$ and $E_t(x_t, y_t)$ are the irradiance distributions of the source and target planes, respectively. Naturally, if the above three constraints can be solved simultaneously, a freeform surface can be obtained to generate the desired irradiance distribution. However, it is very complicated to construct and solve such a set of equations.Instead of directly tackling the above three constraints, the ray-mapping method [3–6], which is the most widely used freeform surface design method, divides the design process into two steps that consider the above three constraints separately. First, the energy conservation equation is solved (Eq. (2)) to obtain the ray-mapping $\mathbf{\phi }: (x_s,y_s)\rightarrow (x_t,y_t)$. Next, by combining the smoothness constraint and Snell’s law, the shape of the freeform surface is constructed point-by-point through numerical integration [3]. However, the ray-mapping that satisfies the energy conservation is not unique, and the main challenge of this method is to obtain an integrable ray-mapping. In the paraxial case, the $L^2$ optimal transport mapping(OT-map) has been well studied for freeform lenses [7–9]. The $L^2$ optimal transport mapping refers to a mapping that satisfies the energy conservation and minimizes the transport cost. The $L^2$ transport cost is defined as follows [10]:
To address the limitations of the OT-map method in non-paraxial and off-axis situations, several modifications have been proposed. For instance, a measure-preserving transformation method is proposed in [15], which alters the initial OT-map iteratively by integrating a Hamiltonian flow. Another approach involves modifying the iterative process of the least-squares solver for OT-map [13,16]. After each iteration, the mapping $m_n$ is used to construct a smooth surface, which is then used to obtain the integrable mapping $m^{'}_{n}$ used in the next iteration. Feng et al. [17] proposed an iterative wavefront tailoring (IWT) method, which considers the smoothness constraint of the wavefront instead of the surface itself and derives a second order partial differential equation for the wavefront. This equation can be iteratively solved to correct the ray-mapping. Importantly, after each iteration, the surface needs to be reconstructed and verified by ray-tracing simulation to determine whether further iteration is necessary. Therefore, these methods essentially iteratively correct the shape of the freeform surface indirectly by correcting the ray-mapping first. Methods that directly correct the shape of a surface based on the results of ray-tracing have not been reported.
In recent years, parameter optimization methods based on differentiable ray-tracing (DiffRT) [18–22] have been extensively studied for spherical and aspherical lenses design in optical imaging systems. In these methods, the lens surface is determined by an analytic expression with a set of variable parameters and its shape can be directly modified according to the results of ray-tracing simulations. However, only a few initial attempts [20,22] have been made on the design of freeform illumination lenses, in which, B-Splines are used to represent the freeform surface and the shape of the surface is optimized directly by updating its control points. These works have chosen a plane or spherical surface as the initial shape, and the design results generally fail to out-perform that obtained with the OT-map method [20].
In this work, in order to take advantages of the OT-map method and overcome its limitations to non-paraxial conditions, we propose an extended ray-mapping method that incorporates differentiable ray-tracing into the design pipeline of the ray-mapping method. The proposed method involves two optimization stages. In the first stage, the freeform surface is preliminarily optimized to match the optimal transport mapping as closely as possible. In the second stage, this surface is further tailored according to the divergence between the target and simulated irradiance distributions directly. Meanwhile, the mean curvature is also constrained in the second stage to make the optimization process more robust and the final freeform surface easier to fabricate. In addition, to enhance the efficiency of the ray-tracing algorithm and increase the design freedom, we use a B-Spline model with uniformly spaced control points and knots, and optimize the calculation process using sparse matrices. Our method enables efficient and accurate design of freeform surfaces that generate complex irradiance distributions under non-paraxial and off-axis situations, as demonstrated through two challenging design examples in this study.
2. Design method
The schematic diagram of a non-imaging optical system is shown in Fig. 1, which generates a specified irradiance distribution with freeform surface. The rays emitted from the source generate an irradiance distribution $I_s(x_s,y_s)$ at the source plane $Z=d_0$. The freeform lens is located at $Z=d_1$ and the rays refracted by the freeform lens generate the irradiance distribution $I_t(x_t,y_t)$ at the target plane $Z=D$. The paraxial condition of the non-imaging optical system means that the maximum deflection angle of rays at the surface is less than $10^{\circ }$ in general. This work focuses on a non-paraxial system and proposes a simple and efficient extended ray-mapping method. In order to implement this method, an analytic expression of the freeform surface and a corresponding efficient differentiable ray-tracing algorithm are necessary.
2.1 B-Spline surface model with uniform control points and knots
B-Spline, which is short for "basis spline", is a mathematical model used to represent arbitrary smooth curves and surfaces, and has been widely used in computer-aided design and 3D reconstruction [23] because of its local controllability. A typical B-Spline surface can be formulated as follows [24]:
where $\theta _{i,j}=(\theta _x, \theta _y, \theta _z)$ is the coordinates of the control points, $m$ is the number of control points along one axis ($m^2$ control points in total), $B_{i,n}(t)$ means the $i-th$ basis function of order $n$, which is a piece-wise polynomial function of degree $n-1$ with a parameter $t$. The basis functions of B-Spline are typically defined and calculated by recursion [24].It is worth noting that it is usually difficult to directly calculate the intersection between the rays and freeform surfaces represented by Eq. (4), which involves obtaining the relationship between the parameters and the real coordinates, i.e., $t_x(x)$ and $t_y(y)$. In order to simplify the solving process, control points and knots are carefully defined, and analytic expressions for efficient ray-tracing are then derived.
Firstly, the knot vector of the B-Spline surface is chosen uniformly, with the interval between adjacent knots as 1. For instance, for a B-Spline surface of order $n$ with $m^2$ control points, the knot vector along one axis is:
Fig. 2. (a) B-Spline surface model with uniformly distributed control points. (b),(c) Visualization results of basis functions and derivative functions of the used B-Spline model of order 3, with 7 control points.
By introducing the above two constraints to the general B-Spline model of order 3, analytic expressions for the basis functions and the derivative functions of the B-Spline surface can be derived. Figs. 2(b) and 2(c) show the visualization results of basis functions and their derivatives with 7 control points. Moreover, with the uniformly spaced control points and knots defined above, an analytic expression for $t$ can be obtained :
It is worth emphasizing that, as shown in Figs. 2(b) and 2(c), for a given value of $t$, only the three adjacent basis functions and their derivatives are non-zero. Therefore, the summation in Eq. (4) actually involves only 9 non-zero basis functions to obtain the surface points. In order to accelerate the simulation, the basis functions and their derivatives are stored as sparse matrices, and the optimized sparse matrix algorithms [25] are used to accelerate the ray-tracing. By utilizing sparse matrices in the simulation, the memory consumption only depends on the number of rays, rather than the number of control points used in the surface model. This enables the use of more control points in the design to provide greater design freedom, without increasing the simulation time and memory consumption.
2.2 Parameter optimization method based on differentiable ray-tracing
With the help of the analytic expression of the B-Spline surface, a parameter optimization method based on differentiable ray-tracing can be implemented. The procedure of the parameter optimization method illustrated in Fig. 3 contains three critical steps.
Firstly, suitable surface parameters $\left \{\theta _i\right \}$ are chosen to establish the structure of the optical system. The positions of the light source and the target plane are determined in accordance with the application requirements. In the second step, rays $\left \{R_i\right \}$ emitted from the source plane are simulated to propagate and refract through the optical system, and finally intersect with the target plane to obtain the intersection points $\left \{P_i(\mathbf{\theta })\right \}$.
Thanks to the analytic expression of the freeform surface, the intersection points between the rays and the freeform surface can be easily calculated using the Newton’s method. The refracted rays at the intersection points can be calculated by first obtaining the analytic expression of the normal vector at the intersection points and then applying Snell’s law. The merit function $m(\mathbf {P})$, which evaluates the performance of the optical system according to the intersection points and is differentiable with respect to the surface parameters, is also calculated in the second step. Notably, the simulation results obtained in this step are consistent with the traditional Monte-Carlo ray-tracing algorithm. The critical difference is that the simulation algorithm in this step is implemented using modern differentiable programming tools such as PyTorch [25]. In the third step, based on the results of differentiable ray-tracing obtained in step two, the gradient between the merit function and the parameters of the surface can be efficiently derived with automatic differentiation [26]. Thereafter, gradient-based non-convex optimizer can be employed to update surface parameters to minimize the merit function.
2.3 Extended ray-mapping method with two optimization stages
Based on the basic theory and implementation details of the optimization method proposed above, the extended ray-mapping method is summarized in Fig. 4.
Fig. 4. Proposed extended ray-mapping method with two optimization stages. In the first stage, the optimal transport mapping is derived and the shape of the freeform surface is tailored to match the ray-mapping according to the deviations between the desired and simulated intersection points. In the second stage, the obtained surface is further tailored according to the difference between the target and simulated irradiance distributions directly.
The proposed method involves two optimization stages. In the first stage, the freeform surface is initialized as a plane and optimized to match the optimal transport mapping as closely as possible. This process is similar to the original OT-map method, except that the traditional point-by-point surface construction method is replaced by a surface construction method based on differentiable ray-tracing, which can obtain the analytical expression of the surface, rather than just the sampling points of the surface. The optimal transport mapping $\mathbf{\phi }$ between the irradiance distributions of the source and target planes is calculated, and two sets of grid points, i.e. $\left \{P^s_i\right \}$ on the source plane (red grids) and $\left \{ P^t_i \right \}$ on the target plane (blue grids) are obtained firstly. Rays $\left \{R_i\right \}$ emitting through $\left \{P^s_i\right \}$ are redirected by the freeform surface and expected to intersect with the target plane at $\left \{ P^t_i \right \}$. These rays are used in the differentiable ray-tracing simulation, and the merit function is defined as the $L^2$ norm between the actual simulated intersection points $\left \{P_i(\mathbf{\theta })\right \}$ and the desired intersection points $\left \{ P^t_i \right \}$ derived by ray-mapping, which can be expressed as follows:
After optimization, the merit function will be minimized and a smooth freeform surface that closely matches the optimal transport mapping will be obtained. Although the freeform surface obtained in the first stage could not generate the desired irradiance distribution accurately in the non-paraxial and off-axis situations, it can be used as a good initial shape for the second stage of optimization.In the second stage, numerous sampled rays and the obtained freeform surface are used in the simulation. The merit function is chosen as the relative root-mean-square difference (RRMSD) between the normalized prescribed and simulated irradiance distributions:
3. Design example
To implement the proposed extended ray-mapping method, we have developed a differential ray-tracing algorithm with B-Spline model using PyTorch [25], an open-source machine learning framework that supports automatic differentiation and GPU acceleration. In addition, in the first stage of optimization, the OT-map is obtained through the back-and-forth method [14], which is a solver for the $L^2$ optimal transport problem with nearly linear time complexity. The following design examples are all executed on a computer with a Ryzen5 5600G CPU and an NVIDIA RTXA4000 GPU with 16 GB of VRAM. By using the B-Spline model introduced in Section 2, our implemented differentiable ray-tracing algorithm costs only about 2.8 s and uses about 6.2 GB of VRAM for one simulation with more than 8 million rays and a $512\times 512$ pixels detector for a B-Spline freeform surface with $512\times 512$ control points.
3.1 Example of non-paraxial illumination design
We first focus on generating complex image pattern under non-paraxial situations. The schematic of the optical system is shown in Fig. 5(a), where the collimated source is located at the plane $Z=0$ mm and emits light of wavelength 650 nm that produces a uniform distribution. The light from this source is refracted by a freeform surface lens and finally intersects the target plane located at $Z=1000$ mm, forming the target irradiance distribution in the region [-466,466] mm $\times$ [-466,466] mm, with a resolution of $512\times 512$, which is generated from a portrait of Carl Friedrich Gauss [29].
Fig. 5. (a) Schematic diagram of the freeform optical system for the first design example. (b) Optimal transport mapping between source and target irradiance distributions.
The freeform lens is a "plane-concave" lens and only the back surface is a freeform surface, which is located at $Z=20$ mm initially. The freeform lens is made of PMMA, with a refractive index of 1.488 at 650 nm, and a dimension of 24 mm $\times$ 24 mm.
In the first stage of the design process, the optimal transport mapping from the light source to the target plane is obtained. The numerical solution size is set as $1001 \times 1001$ and the calculation time is about 16 seconds. The uniform grid points on the source plane, and the corresponding grid points on the target plane obtained through the optimal transport mapping are downsampled and shown in Fig. 5(b).
Next, a set of rays emitted from the uniform grid points on the source plane were used for differentiable ray-tracing simulation. The expected intersection points of these rays with the target plane are the grid points on the target plane mentioned above.
The freeform surface used $512 \times 512$ control points, and a multi-scale optimization strategy is adopted to make the surface optimization process more robust. The Adam [30] optimizer is used for the above optimization process, with learning rates of 0.02, 0.004, and 0.008 used for the three scales($32\times 32$, $128\times 128$ and $512\times 512$), respectively. Fig. 6(a) shows the variation of the merit function during the optimization process. The deviation between the simulated and the target intersection points decreases rapidly, and the shape of the freeform surface remains smooth during the optimization, as shown in Fig. 6(b). Fig. 6(c) shows the target irradiance distribution. The irradiance distortions obtained by ray-tracing simulation after the first optimization stage is shown in Fig. 6(d).
Fig. 6. The first stage of the proposed optimization method for the first design example. (a) The deviation of the target and simulated intersection points decrease during optimization. (b) The shape of the B-Spline surface changed rapidly to minimize the merit function. (c) Target distribution of the first design example (generated from [29]). (d) Simulated irradiance distribution using the surface derived after the first optimization stage. (e) Simulated irradiance distribution using the surface constructed with point-by-point integration method.
It should be noted that in this design example, the maximum divergence angle of the rays exceeds $25^{\circ }$, which does not satisfy the paraxial condition. Therefore, after the first stage of optimization, there is still a certain difference between the simulated and the target irradiance distributions. Fig. 6(e) shows the simulated irradiance distortion using the surface obtained by the traditional point-by-point integration method. It can be seen that due to the large divergence angle of the rays, the irradiance distribution obtained by the point-by-point integration method exhibits obvious distortion, while the freeform surface obtained by the proposed optimization method provides a better match with the target irradiance distribution.
In the second stage of the optimization, 8 million rays are used, and the surface obtained in the first stage is further optimized to minimize the RRMSD between the simulated and target irradiance distributions. As shown in Fig. 7(a), the RRMSD between the simulated and target irradiance distributions gradually decreases from about 20% to 6%, and the simulated irradiance distributions before and after the second stage of optimization are shown in Fig. 7(b). The pixel-wise absolute errors between the target and simulated irradiance distributions before and after the second stage of optimization are obtained. These two pixel-wise errors are normalized together using the maximum error value and visualized in Fig. 7(c). After the first stage of optimization, certain errors around the head and the collar region of the portrait can be discerned, and such errors are significantly suppressed after the second stage of optimization. With a properly chosen initial shape and an efficient differentiable ray-tracing algorithm, the proposed method takes about 8 minutes to obtain the final freeform surface and the second stage of optimization costs only 298 s (2.98 s/it).
Fig. 7. The second stage of the proposed optimization method for the first design example. (a) The RRMSD between the target and simulated irradiance distributions decrease during optimization. (b) Simulated irradiance distributions before (left) and after (right) the second stage of optimization. (c) Normalized pixel-wise absolute errors between the target and simulated irradiance distributions before (left) and after (right) the second stage of optimization.
In addition, we compare the optimization process with and without the use of the curvature penalty term to demonstrate its effect. Fig. 8(a) illustrates the rapid increase in the average absolute value of the mean curvature during the second stage optimization when the curvature penalty term is removed. On the other hand, with the introduction of the curvature penalty term, the mean curvature can be maintained at a lower value consistently. Fig. 8(b) compares the distributions of mean curvature for the optimized freeform surfaces. It is evident that, without the curvature penalty term, the mean curvature exhibits substantial fluctuations, with a maximum absolute value exceeding 3 mm$^{-1}$ (corresponding to a curvature radius less than 0.33 mm). Whereas with the introduction of the curvature penalty term, the mean curvature is effectively controlled within the range of $[-0.19, 0.06]$ mm$^{-1}$, thus ensuring manufacturability.
Fig. 8. The effect of the curvature penalty term. (a) The changes in the average of the absolute value of the mean curvature during the second stage of optimization with and without the curvature constraint. (b) The mean curvature distributions after optimization with (right) and without (left) the curvature constraint.
3.2 Example of off-axis illumination design
In the second example, an off-axis illumination problem with a collimated source is considered. The system setup is shown in Fig. 9(a). A 12-mm-radius collimated source at 650 nm is used, which is located at $Z=0$ mm, with a center point at $(0,0)$ and emitting uniformly distributed rays in the direction of $[0,0,1]$. A "plane-concave" lens with a freeform back surface and coaxial with the source, is located at $Z=20$ mm. The dimensions of the freeform lens is 24 mm $\times$ 24 mm. The target distribution is the images of uniform alphabets "THU" with a resolution of $512\times 512$, which is located at $Z=1000$ mm with the region of [0,400] mm $\times$ [0,400] mm.
Fig. 9. (a) System setup of the second design example. (b) Optimal transport mapping between source and target irradiance distributions. (c) The shape of the B-Spline surface changed rapidly to match the optimal transport mapping. (d) Simulated irradiance distribution using surface constructed by point-by-point integration method. (e) Simulated irradiance distribution using surface constructed with the proposed parameter optimization method.
It should be emphasized that the target irradiance distribution in this example is discontinuous (with regions of zero energy), which makes the optimal transport mapping theoretically contain singular values. It is difficult to achieve high-quality mapping using solvers for continuous optimal transport problem. Constructing surfaces based on such mapping is even more challenging. Fig. 9(b) shows the optimal transport mapping calculated with the back-and-forth method between the uniform source distribution and the target irradiance distribution, and it is evident that the quality of this mapping is not very good. When using the traditional point-by-point integration method to construct the surface, the ray-tracing simulation result shows significant distortion and non-continuous areas, shown in Fig. 9(d). It is due to the combined effect of the off-axis situation and the poor quality of the optimal transport mapping.
Fig. 9(c) shows the optimization process of the freeform surface in the first stage of the proposed parameter optimization method. The shape of the surface gradually changes from the initial plane to the expected tilted surface while maintaining smoothness. Fig. 9(e) shows the result of ray-tracing using the surface obtained after the first optimization stage, which only exhibits weak distortion caused by the off-axis situation. It is evident that the shape obtained through the proposed method is better than that obtained through point-by-point integration.
In the second stage of optimization, 8 million rays are used in ray-tracing simulation to obtain the corresponding output irradiance distribution. The shape of the surface will be directly modified based on the deviation between the simulated and the target irradiance distributions. Fig. 10(a) shows the evolution of the simulated irradiance distributions during the optimization process. It can be seen that the distortion in the output pattern gradually decreases and the irradiance becomes more uniform. The RRMSD between the simulated and target illuminance distributions rapidly decreases from 48% to 12.54%, as shown in Fig. 10(b). Fig. 10(c) further illustrates the distributions of the final simulated irradiance along the line $y=200$ mm on the target plane. Obviously, it is very close to the target value, although there is still some background energy between the letters. Fig. 10(d) shows the final shape of the freeform surface.
Fig. 10. The second stage of the proposed optimization method for the second design example. (a) Simulated irradiance distributions during the second stage of optimization (b) The RRMSD between the target and simulated irradiance distributions decrease during optimization. (c) Target and simulated irradiance distributions along the line y = 200 mm. (d) 3D visualization result and the height map of the desired freeform surface after optimization.
4. Conclusion
In conclusion, an extended ray-mapping method based on differentiable ray-tracing is proposed. By directly adjusting the surface derived by matching optimal transport mapping, this method could extend the OT-map method to non-paraxial and off-axis situations. Additionally, the method combines the advantages of the optimal transport mapping and parameter optimization, achieving a highly efficient and intuitive design process.
In the design examples, only zero-étendue sources are used. However, this method can also be applied to extended sources, and only a ray sampling method for extended sources is required. In addition, more complex optical field control, e.g., controlling both the irradiance distribution and the wavefront at the same time [31,32], can also be implemented by adopting a suitable merit function.
Funding
National Natural Science Foundation of China (61875104, 61927811, 61974080, 61975093, 61991443, 62127814, 62225405, 62235005).
Acknowledgments
The authors thank Zexin Feng for the valuable discussions of the differentiable ray-tracing algorithm and comments on the paper.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. H. Ries and J. Muschaweck, “Tailored freeform optical surfaces,” J. Opt. Soc. Am. A 19(3), 590 (2002). [CrossRef]
2. R. Wu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “A mathematical model of the single freeform surface design for collimated beam shaping,” Opt. Express 21(18), 20974 (2013). [CrossRef]
3. Y. Luo, Z. Feng, Y. Han, and H. Li, “Design of compact and smooth free-form optical system with uniform illuminance for LED source,” Opt. Express 18(9), 9055 (2010). [CrossRef]
4. F. R. Fournier, W. J. Cassarly, and J. P. Rolland, “Freeform reflector design using integrable maps,” Opt. Express 18(5), 5295–5304 (2010). [CrossRef]
5. C. Bösel and H. Gross, “Ray mapping approach for the efficient design of continuous freeform surfaces,” Opt. Express 24(13), 14271 (2016). [CrossRef]
6. R. Wu, Z. Feng, Z. Zheng, R. Liang, P. Benítez, J. C. Mi nano, and F. Duerr, “Design of Freeform Illumination Optics,” Laser Photonics Rev. 12, 1700310 (2018). [CrossRef]
7. X. Mao, S. Xu, X. Hu, and Y. Xie, “Design of a smooth freeform illumination system for a point light source based on polar-type optimal transport mapping,” Appl. Opt. 56(22), 6324 (2017). [CrossRef]
8. K. Desnijder, P. Hanselaer, and Y. Meuret, “Flexible design method for freeform lenses with an arbitrary lens contour,” Opt. Lett. 42(24), 5238 (2017). [CrossRef]
9. Z. Feng, B. D. Froese, and R. Liang, “Freeform illumination optics construction following an optimal transport map,” Appl. Opt. 55(16), 4301 (2016). [CrossRef]
10. C. Villani, Optimal Transport: Old and New (Springer, 2009).
11. 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 (2013). [CrossRef]
12. J.-D. Benamou, B. D. Froese, and A. M. Oberman, “Numerical solution of the optimal transportation problem using the monge–ampère equation,” J. Comput. Phys. 260, 107–126 (2014). [CrossRef]
13. C. R. Prins, R. Beltman, J. H. M. ten Thije Boonkkamp, W. L. 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]
14. M. Jacobs and F. Léger, “A fast approach to optimal transport: the back-and-forth method,” Numer. Math. 146(3), 513–544 (2020). [CrossRef]
15. 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 (2019). [CrossRef]
16. S. Wei, Z. Zhu, Z. Fan, and D. Ma, “Least-squares ray mapping method for freeform illumination optics design,” Opt. Express 28(3), 3811 (2020). [CrossRef]
17. Z. Feng, D. Cheng, and Y. Wang, “Iterative wavefront tailoring to simplify freeform optical design for prescribed irradiance,” Opt. Lett. 44(9), 2274 (2019). [CrossRef]
18. Q. Sun, C. Wang, F. Qiang, D. Xiong, and H. Wolfgang, “End-to-end complex lens design with differentiable ray tracing,” ACM Trans. Graph. 40(4), 1–13 (2021). [CrossRef]
19. Z. Li, Q. Hou, Z. Wang, F. Tan, J. Liu, and W. Zhang, “End-to-end learned single lens design using fast differentiable ray tracing,” Opt. Lett. 46(21), 5453 (2021). [CrossRef]
20. C. Wang, N. Chen, and W. Heidrich, “Do: A differentiable engine for Deep Lens design of computational imaging systems,” IEEE Trans. Comput. Imaging 8, 905–916 (2022). [CrossRef]
21. Y. Nie, J. Zhang, R. Su, and H. Ottevaere, “Freeform optical system design with differentiable three-dimensional ray tracing and unsupervised learning,” Opt. Express 31(5), 7450–7465 (2023). [CrossRef]
22. B. de Koning, A. Heemels, A. Adam, and M. Möller, “Gradient descent-based freeform optics design using algorithmic differentiable non-sequential ray tracing,” arXiv, arXiv:2302.12031v1 (2023). [CrossRef]
23. E. Dimas and D. Briassoulis, “3d geometric modelling based on NURBS: a review,” Adv. Eng. Soft. 30(9-11), 741–751 (1999). [CrossRef]
24. L. Piegl and W. Tiller, The NURBS Book (Springer-Verlag, 1996), 2nd ed.
25. A. Paszke, S. Gross, F. Massa, et al., “Pytorch: An imperative style, high-performance deep learning library,” in Advances in Neural Information Processing Systems 32 (Curran Associates, Inc., 2019), pp. 8024–8035.
26. A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer, “Automatic differentiation in pytorch,” in NIPS-W (2017).
27. C. Wang, “Computational wavefront sensing: Theory, practice, and applications,” Ph.D. dissertation (King Abdullah University of Science and Technology, 2021).
28. J. K. Lawson, C. R. Wolfe, K. R. Manes, J. B. Trenholme, D. M. Aikens, and J. R. Edward English, “Specification of optical components using the power spectral density function,” in Optical Manufacturing and Testing, V. J. Doherty and H. P. Stahl, eds. (SPIE, 1995).
29. C. A. Jensen, “Public domain, via wikimedia commons,” https://commons.wikimedia.org/wiki/File:Carl_Friedrich_Gauss_1840_by_Jensen.jpg (1840).
30. D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv, arXiv:1412.6980v9 (2014). [CrossRef]
31. N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Computation of double freeform optical surfaces using a Monge-Ampère solver: Application to beam shaping,” Opt. Commun. 439, 251–259 (2019). [CrossRef]
32. Z. Feng, D. Cheng, and Y. Wang, “Iterative freeform lens design for optical field control,” Photonics Res. 9(9), 1775 (2021). [CrossRef]
