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

2. The model

We use italics for scalars, boldface for vectors, and $\Vert \cdot \Vert$ 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.

2.1. Givens

We are given an initial base surface described as a height field $z\left(x,y\right)$ 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\left(x,y,z\right)$ which provides the direction $\nabla W$ and divergence ${\nabla }^{2}W$ of the incident rays at any optical surface point $\left(x,y,z\left(x,y\right)\right)$. Slightly abusing terminology, we will refer to ${\nabla }^{2}f$ as the local curvature of the field f. For notational convenience we also define $w\left(x,y\right)$ as the isosurface of W in the neighborhood of $\left(x,y,z\left(x,y\right)\right)$. Finally, we are given a flux density $\overline{s}\left(x,y\right)$ along direction $\nabla W$ at optical surface point $\left(x,y,z\left(x,y\right)\right)$, or equivalently, an incident flux density $s=\overline{s}\cos {\theta }_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 $\nabla 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\left(x,y\right)$ as the pullback of the desired irradiance I and $\widehat{u}\left(x,y\right)$ as the pullback of the actual irradiance $\widehat{I}$ 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 ${\mathbf{r}}_i\in {ℝ}^3$, exitent ray ${\mathbf{r}}_e\in {ℝ}^3$, and refractive index ratio $n\doteq n_i/n_e$ (n = −1 for mirrors), the interface normal n and tangent t satisfy

In our model the surface normal is $\mathbf{n}=\left(-\nabla z,1\right)=\left(-{\partial }_{x}z,-{\partial }_{y}z,1\right)$, the incident direction is ${\mathbf{r}}_i=\nabla W\propto \text{sign}\left(n\right)\left(-\nabla w,1\right)$, and an exit ray rooted at interface point $\left(x,y,z\right)$ meets the projection surface at ray terminus $\left(x,y,z\right)+r{\mathbf{r}}_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 $\Vert {\mathbf{r}}_i\Vert =\Vert {\mathbf{r}}_e\Vert$, in which case Eq. (1) is satisfied by exitent direction

 

Fig. 1 Dilation of a ray bundle from optical surface element dA to projection surface element dB (derived in §2.2–§2.3; see §2.1 for the variable meanings).

Download Full Size | PDF

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: $\int \text{}sdA=\int \text{}udB$. In any area small enough that s and u are essentially constant, we can write the flux constraint as $sdA=udB$. 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\doteq dxdy$ 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, $\frac{1}{2}\left({\mathbf{d}}_1\times {\mathbf{d}}_2\right)$. These vectors are the change in the ray terminus as we move the ray root across the two diagonals of dA: specifically ${\mathbf{d}}_1,{\mathbf{d}}_2\approx dx\left({\mathbf{e}}_x+r{\partial }_x{\mathbf{r}}_e\right)\mp dy\left({\mathbf{e}}_y+r{\partial }_y{\mathbf{r}}_e\right)$ where ${\mathbf{e}}_x$ and ${\mathbf{e}}_y$ 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 ${\mathbf{d}}_1,{\mathbf{d}}_2\approx dx\left({\mathbf{e}}_1+r{\partial }_x\nabla q\right)\mp dy\left({\mathbf{e}}_2+r{\partial }_y\nabla q\right)$, where $\nabla q\doteq \nabla \left(\left(n-1\right)z-nw\right)$ 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 {\theta }_Q=dB\cos {\theta }_B$ where θ_B (resp. θ_Q) is the angle between the local normal of the projection surface (resp. constant-z plane) and the central ray ${\mathbf{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² in Eq. (4) is negligible, so we drop the r² term and divide both sides by $\cos {\theta }_{B}dA$ to obtain the area dilation

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

 

Fig. 2 Algorithmic summary: In the zero-étendue case (left and §3.1–3.2), rays from the incident wavefront propagate through the optic to determine sampling locations on the projection surface. In §3.1, sampled target irradiance intensities are pulled back along the exitent rays to the optical surface and compared with incident intensities from the wavefront to estimate local curvatures for the optic via Eq. (6). In §3.2, both target and actual irradiances are pulled back and compared to obtain a field of curvature corrections via Eq. (8). For positive-étendue light-fields (right and §3.3), the physics are reversed: Pseudo-rays from test points on the target determine sampling locations and directions on a radiant surface. The pushforward of energy flows along theserays forms a system of linear equations in that connects irradiance errors to curvature corrections (Eq. (12)).

Download Full Size | PDF

3.1. 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 ${\nabla }^{2}z$ to reveal the 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 $\left(\overline{s}=o=1,{\nabla }^{2}w=0\right)$, a long projection distance $\left(r\gg 0\right)$, and a nearly flat optical surface $\left(z\approx 0\right)$, 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 $\widehat{u}$ 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,

Algorithmically, we proceed as in §3.1 above, but add calculation of the provisional irradiance $\widehat{u}$ to the simulator’s duties, and update the provisional surface $z\leftarrow 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 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).

 

Fig. 3 Two mirrors estimated via Eq. (8) to produce the same irradiance, but from different base surfaces, as per schematic ray transport diagrams. Grayscale shading on the mirror surfaces indicates local surface curvature.

Download Full Size | PDF

3.3. With light fields, including extended light sources

Thus far we have assumed that the wavefront W associates a single ray vector $\nabla 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 $\mathcal{l}\left(\mu ,\nu ,\varphi ,\theta \right)$ sampled at a (possibly fictive) radiant surface, i.e., $\mathcal{l}$ gives the radiance from surface point $\left(\mu ,\nu \right)$ in direction $\left(\varphi ,\theta \right)$. This contains multiple/extended/asymmetric light sources as special cases; $\mathcal{l}$ 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 $\mathbf{p}\in {ℝ}^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_{\mathbf{p}}dT$, where $I_{\mathbf{p}}$ is the desired average irradiant intensity at p. To wit, we want $I_{\mathbf{p}}dT=\underset{\mathrm{\Omega }}{{\displaystyle \int }}\text{}{s}_{\mathbf{p}}dA$ for the flux density field $s_{\mathbf{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_{\mathbf{p}}$ back through the optical surface z, and calculate a field of dilations $\frac{d{B}_{\mathbf{p}}}{dT}=\frac{dA}{dT}\frac{d{B}_{\mathbf{p}}}{dA}$ 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 ${\mathcal{l}}_{\mathbf{p}}$ of the 4D light field function $\mathcal{l}$ containing those radiance values destined for p. The radiance pushforward field ${\mathcal{l}}_{\mathbf{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:

 

Fig. 4 Axonometric view of an “E”-projecting lens tailored via Eq. (11). The radiant surface (left) simulates two lens-embedded LEDs that subtend between 9° and 22° when viewed from points on the tailored exit surface (middle). The total distance from the radiant to the irradiant surface (right) is 2 lens-widths $=\sqrt{2}\varnothing$. The inset shows the irradiance obtained from a surface designed for a point source, comparably illuminated by the LEDs.

Download Full Size | PDF

We now adapt the dilation-dilution model (Eq. (7)) to find the reverse dilation $d{B}_{\mathbf{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_{\mathbf{p}}\doteq \stackrel{3}{\cos }{\theta }_{Q,\mathbf{p}}\cos {\theta }_{A,\mathbf{p}}$, which echoes the the $\stackrel{4}{\cos }$ law of the camera equation for Lambertian radiant surfaces (for isotropic emitters, ${\mathcal{l}}_{\mathbf{p}}$ will contain a $\stackrel{-1}{\cos }{\theta }_{B,\mathbf{p}}$ term).

Differencing two instances of Eq. (10) — one for a corrected ($c\ne 0$) surface that provides the desired irradiance $I_{\mathbf{p}}$ and one for an uncorrected (c = 0) surface that provides the actual irradiance ${\widehat{I}}_{fp}$ (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 ${\mathcal{l}}_{\mathbf{p}}{o}_{\mathbf{p}}/r_{\mathbf{p}}$ is negligible.

The r.h.s. of Eq. (11) reminds us that the irradiance error is a convolution of the light field restriction ${\mathcal{l}}_{\mathbf{p}}$ and the curvatures ${\nabla }^{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 ${\nabla }^{2}c$ by discretizing the lens surface into facets with area dA and choosing test points ${\mathbf{p}}_i$ on the projection surface. For each test point ${\mathbf{p}}_i$, 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 $N=128\times 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.

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 ${\nabla }^{2}c=f$, can be solved in $O\left(N\log N\right)$ time via FFT. On a rectangular domain, this takes a particularly simple form

$$K\left(c\right)={\mathcal{F}}^{-1}(\mathcal{F}\left(K\left(f\right)\right)\cdot \mathcal{F}\left(\mathcal{L}{)}^{-1}\right).$$

Here $\mathcal{F}(\cdot )$ and ${\mathcal{F}}^{-1}(\cdot )$ are the 2-dimensional DFT and inverse DFT with periodic boundary conditions; $K\left(f\right)$ reflects a field as

$\mathcal{L}$ is the discrete Laplacian operator; $\mathcal{F}{(\mathcal{L})}^{-1}$ is the pseudo-inverse of the Laplacian in the Fourier domain (an elementwise inverse of $\mathcal{F}\left(\mathcal{L}\right)$ with a resulting infinity zeroed out); and ${\mathcal{F}}^{-1}(\cdot )$ is the inverse DFT. All field arithmetic is element wise. The signs in $K(\cdot )$ 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 $1/u\to 1/\left(u+{ϵ}_t\right)$ with small ${ϵ}_t>0$ 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 $\mathbf{Ax}=\mathbf{b}$, the Tikhanov-regularized problem $\left({\mathbf{A}}^{\top }\mathbf{A}+{ϵ}^2\mathbf{I}\right)\mathbf{x}={\mathbf{A}}^{\top }\mathbf{b}$ selects for flatter, smoother corrections.

Appendix B: Variates

Table 1. Variables and their meanings.

View Table

References

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.

3. D. Michaelis, P. Schreiber, and A. Bräuer, “Cartesian oval representation of freeform optics in illumination systems,” Opt. Lett. 36, 918–920 (2011). [CrossRef]   [PubMed]  

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]  

7. R. Wester, G. Müller, A. Völl, M. Berens, J. Stollenwerk, and P. Loosen, “Designing optical free-form surfaces for extended sources,” Opt. Express 22, A552–A560 (2014). [CrossRef]   [PubMed]  

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]  

11. D. Ma, Z. Feng, and R. Liang, “Deconvolution method in designing freeform lens array for structured light illumination,” Appl. Opt. 54, 1114–1117 (2015). [CrossRef]   [PubMed]  

12. M. V. Berry, “Laplacian magic windows,” J. Opt. 1906LT01 (2017). [CrossRef]  

13. W. L. Wolfe, Introduction to radiometry(SPIE, 1998). [CrossRef]