Freeform irradiance tailoring for light ﬁelds

We show how to correct an optical surface to transform an arbitrary incident light ﬁeld into a desired irradiance pattern on a projection surface. Beam dilation errors and optical surface corrections are derived from the pullbacks of the actual and desired irradiances. Étendue eﬀects — the principal obstacle to extended-source tailoring - are factored out by solving a sparse linear system. The method accommodates nontrivial projection surfaces, transport phenomena


Introduction
It is well known that a freeform optical surface can be tailored to convert a zero-étendue beam into almost any desired irradiance field on a projection plane. The problem has attracted diverse modeling and numerical solution strategies, exemplified by, but not limited to, use of PDEs [1], piecewise conic [2] and oval surfaces [3], Monge-Ampère equations [4], shape-from-transport [5], and paraxial caustics [6]. However, modeling assumptions -usually some combination of collimated light, shallow surface slopes, and long effective focal lengths -do not comport with the realities of extended light sources and miniaturized optics. Treatments of the more realistic non-zero étendue problem generally focus on the modest goal of obtaining a uniform irradiance field from an idealized LED, e.g., see [7,8] and all examples in the recent state-of-the-art review by Wu et al. [9]. As noted by Sorgato et al. [10], "The 3D prescribed intensity problem for extended sources with no symmetry has no universal solution yet. " We present a general solution for the irradiance tailoring problem for light fields, a considerably broader class of sources. This enables the design of freeform optics that produce highly structured illumination patterns from high-étendue light sources. As demonstrated in §3.3, the method can produce steeply sloped high-power optics that work at very short effective focal lengths, opening the door to high-efficiency high-brightness systems with very compact dimensions.
Starting with the zero-étendue case, we derive corrections to the optical surface from beam dilation errors, as revealed by comparing the target irradiance with the actual irradiance. In practice, this reduces to a ray-casting and a discrete Fourier transform.
We then generalize to an arbitrary incident light field (plenoptic function). This includes multiple extended light sources of arbitrary shape and directionality. The historical difficulty is that the irradiance is a highly nonstationary convolution of the light field and the optical surface ray-mapping function; this convolution has resisted mathematical characterization that would enable one to invert it and thereby isolate the effects (and thus defects) of the optical surface. Attempts to approximate the deconvolution as a stationary deblurring problem have only shown success in systems with negligible étendue (e.g., [11]). We solve the deconvolution problem by developing a decomposition of the flux through the optic that connects irradiance errors to surface corrections through a system of linear equations. The system is often sparse and rapidly solvable. Where étendue or optical path physics do not allow an exact solution, the method finds a surface that approximates the target irradiance. The resulting method is simple, fast, and versatile -accommodating nontrivial wavefronts, optical paths, and transport phenomena.

The model
We use italics for scalars, boldface for vectors, and · for the Euclidean norm. Most quantities can be interchangeably treated as points, fields, or functions. All variates are tabulated in Appendix 4, and their relationships depicted.

Givens
We are given an initial base surface described as a height field z(x, y) and a target irradiance field I mapped onto a projection surface. Typically z = 0 and I is a bitmapped intensity image. On the radiant side, we assume a wavefront function W(x, y, z) which provides the direction ∇W and divergence ∇ 2 W of the incident rays at any optical surface point (x, y, z(x, y)). Slightly abusing terminology, we will refer to ∇ 2 f as the local curvature of the field f . For notational convenience we also define w(x, y) as the isosurface of W in the neighborhood of (x, y, z(x, y)). Finally, we are given a flux densitys(x, y) along direction ∇W at optical surface point (x, y, z(x, y)), or equivalently, an incident flux density s =s cos θ A , where θ A is the angle of the incident ray with the z-axis. The incident wavefront plays a role in derivations but is cancelled out in the solutions, so it need not be mathematically characterized for the reduction to practice. We assume a transport simulator that follows each ray along ∇W through optical surface z to the projection surface, where it samples the irradiant flux density. A field of these samples, parameterized by where the rays depart the optical surface, is known as the pullback of the irradiance. We will denote u(x, y) as the pullback of the desired irradiance I andû(x, y) as the pullback of the actual irradianceÎ produced by z. For each ray the simulator also provides the z-axis-parallel propagation distance r from the optical surface to the projection surface and the angle of incidence θ B at the projection surface.

Transport from surface gradients
We recall a simple geometric fact governing ray propagation at an interface between two homogeneous isotropic media: Given incident ray r i ∈ R 3 , exitent ray r e ∈ R 3 , and refractive index ratio n n i/n e (n = −1 for mirrors), the interface normal n and tangent t satisfy n ∝ n r i r i − r e r e , equivalently, n t, r i r e = t, r e r i since t, n = 0 .
In our model the surface normal is n = (−∇z, 1) = (−∂ x z, −∂ y z, 1), the incident direction is r i = ∇W ∝ sign(n)(−∇w, 1), and an exit ray rooted at interface point (x, y, z) meets the projection surface at ray terminus (x, y, z) + rr e . Expanding the tangent form of Eq. (1) in these terms and rotating coordinates to the plane of refraction reveals that higher order terms only appear in the Euclidean norms. Dropping those higher order terms is equivalent to assuming that r i = r e , in which case Eq. (1) is satisfied by exitent direction

Dilations from surface Laplacians
We are interested in the bundle of rays that travel from area element dA on the optical surface to area element dB on the projection surface ( Fig. 1). The goal is to tailor the optic such that the actual flux through dA matches the desired flux through dB:´s dA =´u dB. In any area small enough that s and u are essentially constant, we can write the flux constraint as s dA = u dB. Therefore to achieve a desired intensity dilution we must have an equal geometric dilation: This can be reduced to plane geometry in a local Cartesian coordinate frame, where dA dx dy is a square area element with average incident intensity s and dB is a quadrilateral area with average irradiant intensity u: The area of quadrilateral dB is half the cross product of its vector diagonals, . These vectors are the change in the ray terminus as we move the ray root across the two diagonals of dA: specifically d 1 , d 2 ≈ dx(e x +r∂ x r e )∓ dy(e y +r∂ y r e ) where e x and e y are the unit vectors in the xand y-directions. It is more convenient to calculate the area dQ of a quadrilateral formed by taking a constant-z cut through the ray bundle where it meets the projection surface, in which case the diagonals simplify to 2D vectors d 1 , is the exitent ray's rate of lateral displacement as per Eq. (2). For small dA, the two quadrilateral areas are related by dQ cos θ Q = dB cos θ B where θ B (resp. θ Q ) is the angle between the local normal of the projection surface (resp. constant-z plane) and the central ray r e of the bundle. Putting this all together yields the projective area relation If we assume the surface exhibits little curvature inside dA, the coefficient to r 2 in Eq. (4) is negligible, so we drop the r 2 term and divide both sides by cos θ B dA to obtain the area dilation where the second equality uses q = ((n − 1)z − nw) from Eq. (2) and the cosine terms are absorbed into the obliquity term o cos θ B /cos θ Q .

Tailoring
We now develop Eq. (5) into three tailoring fixpoints: In §3.1 we solve directly for the surface z; in §3.2 we eliminate the wavefront w and solve for surface corrections c; in §3.3 we leverage the correction approach into a solution for light-field tailoring. The information flows of these algorithms are summarized in Fig. 2.

As a Poisson problem
Recalling that the intensity dilution (l.h.s. Eq. (3)) must equal the beam dilation (r.h.s. Eq. (5)), we equate them and algebraically isolate ∇ 2 z to reveal the Poisson problem where m so/u − 1 is the mismatch between the light incident on the optic from the source and the light pulled back to the optic from the target. In preparation for future results, we have switched from incident flux density s at dA to directional flux densitys = s/cos θ A , consequently the obliquity term o cos θ A cos θ B sec θ Q gains an additional projective cosine. Treating the Poisson equation (Eq. (6)) as a fixpoint yields a simple tailoring algorithm: (1) Use a provisional surface estimate z to calculate (via simulated transport) the ray-dependent r.h.s. terms in Eq. (6) (target irradiance pullback u, propagation depths r, obliquities o); (2) solve the Poisson problem in Eq. (6) via FFT to update z; (3) repeat until convergence. The fixpoint works well for substantially linear problems where ray bends are modest, and converges quickly if the rays through the initial surface provide a good sampling of the target irradiance pattern. In the essentially linear case of uniform collimated light (s = o = 1, ∇ 2 w = 0), a long projection distance (r 0), and a nearly flat optical surface (z ≈ 0), Eq. (6) contains as special cases freeform methods that view the irradiance as the Laplacian of z, e.g., [6,12] .

Via curvature corrections
A more effective fixpoint is obtained by comparing the desired irradiance u with the irradianceû actually produced by a provisional surface z. We will seek a field c of corrections to surface z that satisfies the linear dilation-dilution model, in which dB c is simply dB in Eq. (5) with c added to z. We write Eq. (7) for u andû (with and without c), difference them, and solve for the Laplacian of c to obtain a curvature correction where dB/dA is the dilation dB/dA in Eq. (5) without the obliquity terms. Note that once we have u andû via simulation, we only need the incident intensities (first form) or dilations (second form) to compute the correction; the wavefront w has been algebraically eliminated. Algorithmically, we proceed as in §3.1 above, but add calculation of the provisional irradiancê u to the simulator's duties, and update the provisional surface z ← z + c. Compute time for each correction scales log-linearly with the number of surface height values, since the dominant calculation is the DFT in the Poisson solver. We typically start with the simplest surface that provides good coverage of the target; convergence is fastest when most of the rays sample most of the target density. For example, tailoring a mirror represented by 128 × 64 height values to accurately reproduce a detailed high-contrast image of a human eye (Fig. 3) took less than one second on a laptop computer. See Appendix 4 for numerical considerations.
The correction method has some advantages: Most of the nonlinearities of the optical path are captured on the r.h.s. of Eq. (8) and can be calculated exactly in the simulation of transport, while the linearizations of §2.2- §2.3 are at play only in the relatively small curvature corrections, where they are accurate. Our experience is that Eq. (8) converges faster than Eq. (6) for "easy" problems and yields superior results for harder problems where steep slopes, high curvatures, short effective focal lengths, Fresnel losses, or high contrast ratios are present. An added utility is that the tailored optic can inherit some desired properties from a suitably designed base surface, e.g., demagnifying to increase brightness or crossing all rays to obtain globally convex or concave surfaces (e.g., the r.h.s. of Fig. 3).

With light fields, including extended light sources
Thus far we have assumed that the wavefront W associates a single ray vector ∇W to every point in space. While more general than a collimated or point light source, this falls far short of nature's full plenoptic function. For the purposes of tailoring, the plenoptic function can be summarized as a light field of spatially-and directionally-varying radiances (µ, ν, φ, θ) sampled at a (possibly fictive) radiant surface, i.e., gives the radiance from surface point (µ, ν) in direction (φ, θ). This contains multiple/extended/asymmetric light sources as special cases; can also incorporate the effects of other surfaces prior to the freeform in the optical path. We note that the second law of thermodynamics almost always excludes the possibility of an exact solution for a tailoring problem with a positive-étendue source, but good approximations have practical utility.
To extend the correction method to incident light fields, we consider both the radiant and irradiant surfaces, but reverse their roles. We decompose the flux through the optic into a superposition of "light-collecting" spherical wavefronts that propagate backward from the projection surface. Adapting the dilation-dilution model (Eq. (7)) to relate light collection in each wavefront to the freeform's local curvatures, the superposition becomes a sparse system of linear equations whose solution isolates and quantifies the freeform's curvature errors.
To work out the details, it is useful to imagine a camera that views the optical surface through a small aperture of area dT located at a point p ∈ R 3 on the projection surface. The camera sees a distorted image of the light source(s) on a subset Ω of the optical surface. Ideally, the total observed flux matches I p dT, where I p is the desired average irradiant intensity at p. To wit, we want I p dT =´Ωs p dA for the flux density field s p that the camera at p observes on the optical surface. We can calculate this quantity from the area dilation (Eq. (5)) by reversing the roles of the radiant and irradiant surfaces: We imagine projection surface point p emitting a spherical wavefront w p back through the optical surface z, and calculate a field of dilations d A of this wavefront as it propagates backward from the projection surface to the optical surface z and then to the radiant surface, where we take a 2D subset p of the 4D light field function containing those radiance values destined for p. The radiance pushforward field p plays an analogous role to the irradiance pullback field u in previous sections and is similarly parameterized by the optical surface's coordinate system. Conservation of energy requires that each step have the same total flux: where the integrand on the right-hand side is the standard radiometric calculation of flux from a radiant patch dB p through an aperture dA, for each dB p and dA pair that deliver light to dT [13]. We now adapt the dilation-dilution model (Eq. (7)) to find the reverse dilation dB p . First, we are now tracing light backward, so we must negate z + c and invert n. Second, the dilation of interest is now between the optical surface to the radiant surface; the relationship between the optical surface "facet" size dA in Eq. (7) and the projection surface "pixel" size dT is arbitrary. For simplicity of exposition we choose dA = dT. Applying these changes to Eq. (7), solving it for dB, and substituting the result into the irradiance integral (Eq. (9)) yields a relation between curvature corrections c and the irradiance at p: Here all projective cosines have been absorbed into the obliquity term o p cos 3 θ Q,p cos θ A,p , which echoes the cos 4 law of the camera equation for Lambertian radiant surfaces (for isotropic emitters, p will contain a cos −1 θ B,p term). Differencing two instances of Eq. (10) -one for a corrected (c 0) surface that provides the desired irradiance I p and one for an uncorrected (c = 0) surface that provides the actual irradianceÎ p (as per the r.h.s. of Eq. (9)) -yields the irradiance error: Here we have assumed that the correction is small enough that the change to p o p /r p is negligible. The r.h.s. of Eq. (11) reminds us that the irradiance error is a convolution of the light field restriction p and the curvatures ∇ 2 c of the correction field. That convolution is nonstationary because the p-subscripted terms are different for each test point. However, unlike the irradiance convolution, the curvature-correcting deconvolution is straightforward: We set up a least-squares problem for the Laplacian of the correction ∇ 2 c by discretizing the lens surface into facets with area dA and choosing test points p i on the projection surface. For each test point p i , Eq. (11) where the first vectorized field is nonzero wherever the camera at p sees light coming from the optical surface z. Iterating over all test points p i , this forms a large but sparse system of linear equations for the elements of ∇ 2 c. To tailor the surface z we solve the system of linear equalities for ∇ 2 c, solve the Poisson problem for the correction c, add c to surface z, and repeat until convergence. Figure 4 shows a lens with two embedded LEDs that was tailored in this manner to project a white "E" on a black background. For N = 128 × 128 height values the calculation took 25 seconds, and can be expected to grow at a rate between log-linearly (solution time is dominated by the Poisson problem) and cubically with N (solution time is dominated by the linear system), depending on how much of the optical surface is involved in irradiating each test point.
In general, the positive-étendue tailoring problem is a strongly non-convex optimization problem with no zero-error solution and many local minima. Good results depend on favorable initializations. We often find that a zero-étendue design is a useful starting point for high-étendue problems; such is the case in Fig. 4, where a high-fidelity design for a point light source produces a highly blurred irradiance but serves as a good initial surface for high-étendue tailoring.

Summary and open questions
Observing an essentially linear relationship between incident wavefront curvature, optical surface curvature, beam dilation, and irradiant intensity dilution, we formulated zero-étendue freeform irradiance tailoring as a Poisson problem. Algebraically eliminating the wavefront yielded a sequence of surface curvature corrections based on pullbacks of the irradiance; these rely on considerably more modest linearizations than the basic Poisson problem, and perform well in optical path geometries previously considered difficult. Reversing the roles of the radiant and irradiant surfaces in this formulation, we found a solution for general positive-étendue light-field tailoring, e.g., with multiple extended sources. In particular, the problem of discovering what parts of an optical surface are responsible for flaws in the irradiance -previously a poorly understood nonstationary deconvolution -is revealed to be a straightforward system of linear equations in the domain of curvature correction fields. This approach easily accommodates varied light transport phenomena; beyond the pedagogical examples in this paper, we have used it to tailor freeforms situated in optical paths with spatially separated extended light sources, additional optical elements, partial occlusions, internal reflections, and Fresnel losses.
Strictly speaking, all realistic tailoring problems involve positive-étendue sources, and almost all positive-étendue tailoring problems are infeasible -existence of exact solutions would imply entropy-reduction. Yet we find that many problems have good approximate solutions. Since approximation quality can be improved by discarding some radiance to reduce étendue, an open question of great interest is how to characterize and optimize this trade-off for specific irradiance targets and optical path geometries.