New approach to construct freeform surface by numerically differential formulation

Abstract. We proposed a new design method for single freeform reflective (or refractive) surface tailored to redistribute the radiant flux onto a prescribed illumination pattern. Unlike the conventional optimization approaches based on the grid mapping, in this study we estimated each segmental freeform surface by locally solving a second-order differential equation, which formulates the energy transportation between each domain cell. With finite element method via Hermite element, we validated a series of smooth reflective/refractive surfaces to reallocate the radiant flux from a point source toward a target plane with specific patterns. The proposed technique offers a large flexibility by varying the vectors of each ray with multiple refraction (or reflection), which imposes no restriction on the target distribution, collective solid angle, or even target topography.


Introduction
Construction of a freeform surface to reallocate the radiant flux is a substantial challenge in terms of nonimaging optics. The basic concept is to convert the radiant flux emitted from a point source onto a prescribed illumination, usually in Cartesian coordinate. Parkyn 1 firstly pioneered a tessellation method to map the cell area between a point source in polar grid with a particular rectangular illumination in Cartesian grid. The cells of both grid topologies were sized so that they share common energy. The methodology was based on the assumption that when the polar cell on the light source was exactly mapped to the corresponding rectangular target cell, the flux of the angular cell would be fully transferred to the target cell accordingly. Each corresponding pair of cells on these grids has an associated surface normal vector, and then governs the surface shape using extrinsic differential geometry. The major concern about equal flux grid method lies in a fact that the constructed surface, either reflective or refractive, is unlikely to guarantee the integrability condition in the overall surface, and thus leads to a discontinuous boundary. The cumulative deviations from the ideal surface normal become more pronounced at the far corners of the output grid. The discontinuous boundary could lead to a problem. The problem is that a small amount of radiant flux will be distorted by the sharp cliff. Hence, sufficient partition cells shall be numbered to confirm the estimated accuracy for surface construction.
To alleviate the discontinuous issue of freeform, Wang et al. 2 introduced the variable separation mapping strategy to take into account the normal deviation of the curve below a present value. The surface curve was regenerated with tangent vectors along the curve perpendicular to the normal vector, thus eliminating the accumulated normal deviation. The proposed scheme constructs a smooth surface by interpolating the discontinuity surface with non-uniform rational basis spline (NURBS). 3 On the other hand, the integrability condition η ¼ N · ð∇ × NÞ ¼ 0 is necessary in order to make a smooth freeform surface, where N is the surface normal. 4 To fulfill the integrability condition, Fournier 3 interpolated the ellipsoids, obtained from Oliker's algorithm with appropriate grid partition, by using the NURBS where parameters of the NURBS were determined by minimizing the residual curl η. Although the surface can be made smooth and the resulted radiant flux can be less distorted via the above mentioned approaches, the amount of energy deviation is still difficult to be estimated before the freeform surface is constructed. Furthermore, the computational costs of the constrained optimization for the NURBS' parameters may rise significantly and the convergence is not guaranteed as the number of grid partition increases.
In additional to the tessellation method, the freeform surface design problem is also related to the solution of the wellknown Monge-Ampere (MA) equation. In 1972, Schruben first derived the MA equation obtained from reflector design. Beginning in 1980s until the present, Oliker et al. formulated the freeform problem into a highly nonlinear MA type of partial differential equation (PDE). [5][6][7][8][9] The MA PDE also appears in other research fields, such as differential geometry, optimal control and mass transportation, meteorology, and geostrophic fluid. [10][11][12] Although many works about existence and regularity of MA PDE have been conducted, numerical computation of freeform surface reconstruction by solving the MA PDE remains a challenge in practice. We are curious if there exists a systematic way to construct the freeform in which not only the surface is guaranteed to be smooth, but also the correspondence between light source and target plane is ensured under a controllable energy deviation.
In this article, we formulate the ray forward and inverse propagation with negligible energy deviation within the MA PDE framework. By following the works of Caffarelli and Oliker 13 , the freeform considered here is assumed to be convex. The model is based on local energy conservation in which the flux correspondence associated with each local domain is prescribed. We take advantage of mathematical analysis of the MA PDE and successfully formulate the freeform design problem into the popular finite element paradigm. The formulation could maintain the ability to control the light beam and guarantee a smooth freeform surface. Moreover, this article describes this in detail by going through the different modeling steps of the model. All the main restrictions influencing the target distribution, collective solid angle, or source distribution are discussed, alongside the viability of realizing an optical surface for different specific requests.

Freeform Model in Localized PDE
To develop a method for freeform construction subject to the source-target information, we must deduce an adequate, physically realistic model of ray propagation in the free space. The model can then be employed to understand the reflective (refractive) forward process of rendering illuminance according to a known optical configuration. Also, the model is applicable to the inverse algorithm of estimating the freeform surface from a given source-target corresponding condition. First of all, the ray transportation must satisfy the energy conservation among coordinate transformation. 14,15 The radiant intensity I (Watts∕sr) of a point source (in polar domain Ω θ ) and prescribed illuminance distribution E (Watts∕m 2 ) of the target region (in Cartesian domain Ω x ), in one dimension, can be related as Z The goal is to design a freeform surface by analyzing the potential and limitation of using the PDE. In this article, the freeform surface is assumed to be continuous and convex. By partitioning the global domain Ω θ and Ω x into many local subdomains Ω i θ and Ω i x , so that each corresponding pair satisfies the equal flux i ¼ 1; 2; : : : ; n, we are able to ensure the convergence of the integration. As long as the number of partitions is sufficiently large, the freeform can be considered to be a flat surface. As a result, a reduced but simplified representation of a freeform can be derived on each local domain. Hereafter, a flat surface criterion to the freeform in each local domain is called a small local planar approximation.
The representation of the local freeform is indeed a nonlinear differential equation and shall be derived in the following. First, let us introduce some notations. Global coordinate frame is denoted by ½ 0 with basis vectors fe 0 1 ; e 0 2 g and origin point O. The target plane T is expressed by a plane equation and ½X 0 are the normal vector, shift and spatial variables, respectively.
Although the laws of reflection (or refraction) are invariant with coordinate system, solution of freeform problem requires that the relations derived from these laws can be expressed in a coordinate system appropriate to the geometry of the given problem. Hence, we define a local coordinate system, by which the i'th local coordinate of any vector V is denoted by ½V i with corresponding orthogonal basis fe i 1 ; e i 2 g and local origin O i . Coordinate transformation from i'th to j'th coordinate system can be conducted via a sequence of operators, where R j i and T j represent the rotation of basis vectors and translation: respectively. The scenario of our freeform construction is schematically in Fig. 1. The segmental freeform surface to be constructed on i'th local coordinate system can be analytically presented as a function ½A i ¼ ½u; fðuÞ with surface normal as As shown in Fig. 1(a), while a ray from point source O is hitting on the point A at the freeform reflector, the ray is guaranteed to be transported toward the target plane T. The refraction law in vector form is where n 1 and n 2 represent the refraction indexes near the source side and target plane side, respectively. For the reflection law, 7 the refraction index can be defined by n 1 ¼ 1.0 and n 2 ¼ −1.0. Here, we size each cell of subdomain pair (Ω i θ and Ω i x ) to allow equal flux transfer. In order to derive the freeform surface via the preceding equations, we must introduce an intermediate variable u in the computational domain (Ω i u ) to link the computation in both source and target domains (Ω i θ and Ω i x ). By changing the variable, Eq. (3) can be rewritten as Obviously, the integral form of energy conservation, Eq. (1) can be represented by a differential equation in summation: 7 Let ðu i h ; 0Þ be the coordinate of the hitting point where the incident ray OA ! with light source ½O i ¼ ðL i x ; L i z Þ hits on i'th local x-axis. The x-coordinate of the hitting point u i h can be expressed as two separate functions of u and θ i by trigonometric formula: where θ i and ϕ i 0 are the angle of OO !i with OA ! and e i 2 , respectively, and ρ 0 is the distance between O and O i . The transformation between variables u and θ i can be linked by the joint point u h and the Jacobian chain rule becomes Similarly, suppose the incident ray OA ! is reflected onto the point ½X i ¼ ðx i ; z i Þ at the target plane T. The variables x i and u related by the law of reflection can be expressed as Dividing Eqs. (9) by (10), the equation that defines a local freeform surface can be obtained: If we fix the source and target plane position, the freeform surface defined in Eq. (11) can also be considered as a parameterized surface where the parameter is the x-coordinate of the i'th coordinate system. The freeform parameters depend on the external geometry at which each optical route O-A-X constitutes an epipolar plane. In order to transform the underdetermined Eq. (11) to a determined problem, we replace z i by x i through the equation C 0 · ℋ 0 i ½X i ¼ d 0 : where Rðu; fÞ ¼ Oðu 2 ; uf; uf 0 ; f 2 ; ff 0 ; f 02 Þ is the high order nonlinear term that can be neglected under the small local planar approximation, and the coefficients, C i j , j ¼ 1; 2; : : : ; 7, are constants with respect to the i'th local coordinate system. For the refraction case, the variables x i and u can be related by the law of refraction following a similar derivation. After differentiating Eq. (12) to obtain Jacobian dx i ∕du and inserting Eqs. (8) and (12) into Eq. (5), we can devise the local reflection (refraction) freeform surface to a nonlinear differential equation: where α i 1 , α i 2 , α i 3 , and α i 4 are coefficient functions subject to low order derivatives of f. This equation describes the raytracing map and the energy redistribution. The total energy is conserved as ∫ Ω i x EðuÞxðuÞdu in each local domain. Therefore, the energy deviation is mainly contributed by the neglected term Rðu; fÞ and can be controlled when the integral ∫ Ω i x EðuÞRðu; fÞdu can be estimated. The construction of freeform surface by numerically differential formulation has an advantage in that we can regularize the smoothness tolerance by setting the partition number (Ω i θ and Ω i x ) in a tradeoff between computational efficiency and surface smoothness. In order to ensure the numerically estimated accuracy of freeform surface, sufficient partition cells shall be numbered in both domains. Here, we impose the edge rays R i−1 and R i in the polar cell toward the corresponding boundary ½b i−1 ; b i in the Cartesian cell, so the radiant flux is fully transferred without stray loss. The edge ray,R i , reflected by the freeform surface is propagating toward the corresponding point b i . We shall set a reasonable right end point of the intermediate cell Ω i u , by which the freeform possesses convexity in the local domain. To solve the nonlinear PDE in Eq. (13), we employ the vanished moment method proposed by Feng and Neilan and use Newton's method to find the adequate root function. The linearized equation at each iteration step can be solved by the finite element method. Here, we skip some straightforward but tedious expansion and rearrangement. Figure 2 demonstrates the key steps of the freeform construction process. The algorithm of freeform construction in this session can be summarized as following:

Freeform Surface via Local Reconstruction
Algorithm: Error Evaluation In order to validate the proposed freeform scheme, we first constructed a surface by the previous algorithm and compared it with a known parabolic reflector; both are aimed to deliver a collimation beam as shown in Fig. 3. A point source was placed at the focal point ðx; yÞ ¼ ð0; 10Þ of an ideal parabolic reflector y ¼ −x 2 ∕12 þ 13. Here, we limit the emitted angle from the point source within the range of Ω θ ¼ ½−119 deg; 119 deg, and the radiant flux reflected by the ideal parabolic reflector was directed toward the given target plane Ω x ¼ ½−10; 10 with a uniform and parallel propagation. Figure 3(a) exhibits the surface comparison between the curve constructed via proposed algorithm and ideal parabolic curve. The surface error e s ¼ max x∈Ω x kS FF ðxÞ − S P ðxÞk, maximum deviation between the constructed freeform surface (S FF ) and the ideal parabolic curve (S P ), is on the order of 10 −4 to 10 −6 , which is extremely small and quadratically inverse to the number of partition cells. The results are in agreement with the theoretical prediction that the error comes from the finite element method with Hermite element. [16][17][18] 4 Construction Algorithm: Two Design Types

Uniform Illuminance
After validating our freeform algorithm with a comparative parabolic surface by forward testing, we now aim to achieve a uniform illuminance with a geometric correspondence by the freeform surface via local reconstruction algorithm (FS-LRA). As shown in Fig. 4, a 20-m-wide target plane was placed beneath a point source with 6-m distance. The emitted angle from the point source was confined for reflection and refraction within Ω θ ¼ ½−30 deg; 30 deg. With a mere 80 partition cells and 960 partition cells in both source-target subdomains, a freeform reflector and refractor can be generated in about 1.1 and 10 s, with 2.3 GHz Intel® dual core computer. The study employed a ray tracing simulation on the constructed freeform surface with a commercial tool LightTools®, in which light rays from the source were incidental to the target plane. The root-mean-square (RMS) nonuniformity is merely 1.05% (reflector) and 2.10% (refractor), respectively. It is noted that the accuracy or partition number of FS-LRA is highly dependent on the tolerance of the freeform surface. For uniform illuminance request in both Step 2: The global information including the location of light source and the target plane is transformed into i'th coordinate.
Step 3: The computational domain Ω i u is determined by the corresponding subdomain pair ðΩ i θ ; Ω i x Þ, the refraction or reflection laws and the convexity of the freeform.
Step 4: The freeform surface is calculated by solving Eq. (13) on Ω i u using Hermite finite element method.
Step 5: The (i þ 1)'th coordinate system is defined by the rightmost boundary point O iþ1 and the normal and tangent vectors at O iþ1 . If i is not equal to n, go to step 2.
Step 6: The freeform surface is assembled from local surfaces. reflect/refractive cases, a monotonic curve is accessible with the prescribed simple target surface. Because the mathematical form in the refraction case is more complex than that in reflection case, 960 partitions are required for convergence to achieve the design uniform illumination pattern, but only 80 partitions are required for the constructed reflector to achieve the same uniformity.

Linear Varying Illuminance
In this example, we extended our technique to create a prescribed color illumination by a pair of independent sources with different colors. Figure 5 shows the schematic diagrams of the system layout for both reflect/refractive cases. With a 20-m target plane for both cases, we aimed to create a linearly varying international commission on illumination (CIE) chromatic ðx; yÞ along the x-direction associated with two color sources. The distance between source and target plane is 10 and 70 m, respectively. The light sources leave the center at 10 m and 15 m, respectively. The different configuration conditions of reflection and refraction design are to achieve desired illuminance performance. As shown in Figs. 5(c) and 5(d), the illuminances feature a high linearity along the x-direction. There are only about 5.58% and 1.52% RMS nonuniformity with the perfect chromatic line for reflection and refraction cases, respectively.

Conclusion
In this article, we proposed a freeform model based on local energy conservation associated with a series of coordinate transformation. The major advantage of this method is to exploit the well-proven existence and continuity of an elliptic nonlinear differential equation of the Monge-Ampere type. We successfully formulated the freeform design into a finite element paradigm with Newton's iteration. The key maneuver is to unify the variables (θ and x) in both domains (Ω i θ and Ω i x ) into a common intermediate variable u in the computational domain (Ω i u ). After imposing the small local planar approximation, which is absolutely valid with a sufficient number of partition cells, a complete freeform reconstruction algorithm can be developed. In addition, the segmental differential formulation is capable of taking numerical error into account in each step. The proposed surface parameterization guarantees the surface smoothness with no restriction on the desired target distribution, collective solid angle, or source distribution, respectively. This technique still leaves many opportunities open and clearly more research must be carried The maximum surface difference between ideal parabolic curve and numerical surface. The surface error with respect to subdomain numerical approximates a function Sðx Þ≈ 1∕x 2 . The amount of error reduces quadratically as we increase the number of partition cells. Fig. 4 (a) and (b) Layout of the freeform surface reflector and refractor is designed for uniform illuminance, respectively. (c) Illuminance distribution on the target plane. Inset is the magnification of illuminance plateau pattern. The nonuniformity of reflection case was well controlled within 1.05% RMS and 1.89% peak-to-valley. (d) The nonuniformity of refraction case was well controlled within 2.10% RMS and 4.01% peak-to-valley. out to explore its potential in full. First, at this moment, the issue was tackled in a one-dimensional case; complete treatment in consideration of twist deformation in freeform surface is underway. Second, point source approximation is another subject to be addressed when the freeform structure was placed in proximity of the light source such as light emitting diodes applications. The preliminary results presented in this article, however, indicate that parameterizing the freeform surface via local freeform PDE may create another route in the field of freeform optics, as it has the potential to take freeform design into many nonimaging applications.