We show how to correct an optical surface to transform an arbitrary incident light field 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 effects — the principal obstacle to extended-source tailoring — are factored out by solving a sparse linear system. The method accommodates nontrivial projection surfaces, transport phenomena, and incident wavefronts, including those from multiple extended light sources. Numerical experiments achieve high fidelity and contrast ratios in as little as -time for a surface represented by N height values.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
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 , piecewise conic  and oval surfaces , Monge-Ampère equations , shape-from-transport , and paraxial caustics . 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. . As noted by Sorgato et al. , “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., ). 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.
2. 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 B, Table 1 and their relationships depicted in Fig. 1.
We are given an initial base surface described as a height field 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 which provides the direction and divergence of the incident rays at any optical surface point . Slightly abusing terminology, we will refer to as the local curvature of the field f. For notational convenience we also define as the isosurface of W in the neighborhood of . Finally, we are given a flux density along direction at optical surface point , or equivalently, an incident flux density , 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 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 as the pullback of the desired irradiance I and 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.
2.2. Transport from surface gradients
We recall a simple geometric fact governing ray propagation at an interface between two homogeneous isotropic media: Given incident ray , exitent ray , and refractive index ratio (n = −1 for mirrors), the interface normal n and tangent t satisfy
In our model the surface normal is , the incident direction is , and an exit ray rooted at interface point meets the projection surface at ray terminus . 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 , in which case Eq. (1) is satisfied by exitent direction
2.3. 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: . In any area small enough that s and u are essentially constant, we can write the flux constraint as . 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 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 where and are the unit vectors in the x- and 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 , where is the exitent ray’s rate of lateral displacement as per Eq. (2). For small dA, the two quadrilateral areas are related by where θB (resp. θQ) is the angle between the local normal of the projection surface (resp. constant-z plane) and the central ray 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 r2 in Eq. (4) is negligible, so we drop the r2 term and divide both sides by to obtain the area dilationEq. (2) and the cosine terms are absorbed into the obliquity term .
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.
3.1. As a Poisson problem
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 , a long projection distance , and a nearly flat optical surface , Eq. (6) contains as special cases freeform methods that view the irradiance as the Laplacian of z, e.g., [6, 12].
3.2. 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,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 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 irradiance to the simulator’s duties, and update the provisional surface . 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 A 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).
3.3. With light fields, including extended light sources
Thus far we have assumed that the wavefront W associates a single ray vector 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 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 , where is the desired average irradiant intensity at p. To wit, we want for the flux density field 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 back through the optical surface z, and calculate a field of dilations 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 of the 4D light field function containing those radiance values destined for p. The radiance pushforward field 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:13].
We now adapt the dilation-dilution model (Eq. (7)) to find the reverse dilation . First, we are now tracing light backward, so we must negate 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 . 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 , which echoes the the law of the camera equation for Lambertian radiant surfaces (for isotropic emitters, will contain a term).
Differencing two instances of Eq. (10) — one for a corrected () surface that provides the desired irradiance and one for an uncorrected (c = 0) surface that provides the actual irradiance (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 is negligible.
The r.h.s. of Eq. (11) reminds us that the irradiance error is a convolution of the light field restriction and the curvatures 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 by discretizing the lens surface into facets with area dA and choosing test points on the projection surface. For each test point , Eq. (11) generates one linear equality constraint, in inner product form
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 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.
4. 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.
Appendix A: Numerical considerations
The Poisson problem , can be solved in time via FFT. On a rectangular domain, this takes a particularly simple form
Here and are the 2-dimensional DFT and inverse DFT with periodic boundary conditions; reflects a field as
is the discrete Laplacian operator; is the pseudo-inverse of the Laplacian in the Fourier domain (an elementwise inverse of with a resulting infinity zeroed out); and is the inverse DFT. All field arithmetic is element wise. The signs in are optionally positive to impose Neuman boundary conditions on the corrections (zero derivative in the direction normal to the boundary); this fixates the boundary of the projection in the case of a collimated source. Tailoring can also be done in polar or spherical coordinates using mixed Fourier/Chebyshev methods for the Poisson problem.
Pulled-back zeroes and near-zeroes in the target irradiance field can imply infinite local dilation and thus infinite curvature. To prevent numerical instabilities, we recommend the substitution with small decreasing in each iteration of the fixpoint.
The linear system in light field tailoring Eq. (11) can have poor numerical condition or be underconstrained, e.g., if the irradiance is undersampled or if some optical surface area element is not on any path from a lightsource to a projection test point. In such cases, instead of , the Tikhanov-regularized problem selects for flatter, smoother corrections.
Appendix B: Variates
1. H. Ries and J. Muschaweck, “Tailored freeform optical surfaces,” J. Opt. Soc. Am. A 19, 590–595 (2002). [CrossRef]
2. V. Oliker, Trends in nonlinear analysis(Springer, 2003), chap. 3.
4. C. Prins, J. ten Thije Boonkkamp, J. van Roosmalen, W. IJzerman, and T. Tukker, “A numerical method for the design of free-form reflectors for lighting applications,” Tech. Rep. CASA-report 1322 (2013).
5. Y. Schwartzburg, R. Testuz, A. Tagliasacchi, and M. Pauly, “High-contrast computational caustic design,” ACM Trans. Graph. 33, 74 (2014). [CrossRef]
6. G. Damberg, J. Gregson, and W. Heidrich, “High brightness hdr projection using dynamic freeform lensing,” ACM Trans. Graph. 35, 24 (2016). [CrossRef]
8. Z. Zhao, H. Zhang, H. Zheng, and S. Liu, “New reversing freeform lens design method for led uniform illumination with extended source and near field,” Opt. Commun. 410, 123–129 (2018). [CrossRef]
9. 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, 1700310 (2018). [CrossRef]
10. S. Sorgato, R. Mohedano, J. Chaves, M. Hernández, J. Blen, D. Grabovičkić, P. Benítez, J. C. Miñano, H. Thienpont, and F. Duerr, “Compact illumination optic with three freeform surfaces for improved beam control,” Opt. Express 25, 29627–29641 (2017). [CrossRef] [PubMed]
12. M. V. Berry, “Laplacian magic windows,” J. Opt. 1906LT01 (2017). [CrossRef]
13. W. L. Wolfe, Introduction to radiometry(SPIE, 1998). [CrossRef]