Freeform Optics: Optimal Transport, Minkowski Method, and Monge–Ampère-Type Equations

1. Description and Background

A freeform optical surface, simply stated, refers to a surface whose shape lacks rotational symmetry. The use of such surfaces allows generation of complex, compact, and highly efficient imaging systems. Ever since lenses without symmetries were used in World War I in periscopes, the engineering and design of freeform optical surfaces have gone through remarkable evolution, with applications in a wide range of areas, including medical devices, clean energy technology, military surveillance equipments, mobile displays, remote sensing, and several other areas of imaging and nonimaging optics that can benefit from distributing light from a source to a target in a controlled fashion using spatially and energy-efficient systems. See 17 and the references therein.

Mathematically, the design of freeform optical surfaces is an inverse problem related to optimal transportation theory and leads to a class of nonlinear partial differential equations (PDE) called generated Jacobian equations for which the Monge–Ampère equation is a prototype. The modeling of the problem is based on the systematic application of the laws of reflection and refraction in geometric optics along with energy conservation principles.

For completeness we recap the laws of geometric optics. Consider a reflecting or refracting surface $\mathcal{S}$ and suppose a unit vector in the direction of $x$ is incident on the surface $\mathcal{S}$ at a point where the normal to $\mathcal{S}$ is given by the unit vector $\nu$. Reflection and refraction happen in tandem. (Fig. 1). However, for simplicity of the models we treat them separately. If $r$ is a unit vector in the direction of the reflected ray, the law of reflection (angle of incidence is equal to angle of reflection) in vector form says:

$$\begin{equation} r = x - 2\langle x, \nu \rangle \nu . {\tag{1.1}} \end{equation}$$

If medium $I$ and medium $II$ are two homogeneous isoptropic media with respective refractive indices $n_1$ and $n_2$ and $m$ is the direction of the refracted ray into medium $II$, the law of refraction (Snell’s law), which states $\dfrac{\sin \theta _i}{\sin \theta _t} = \dfrac{n_2}{n_1}$ where $\theta _i$ is the angle between $x$ and $\nu$ (angle of incidence) and $\theta _t$ is the angle between $m$ and $\nu$ (angle of refraction) in vector form, states that:

$$\begin{equation} m = \dfrac{1}{\kappa }(x - \lambda \nu ) {\tag{1.2}} \end{equation}$$

where $\kappa = n_2/n_1$ is the relative refractive index, $\lambda = \lambda (\kappa , x, \nu )$ and $\nu$ points into medium $II$.

Figure 1.

Laws of reflection and refraction.

When medium I is optically denser than medium II or equivalently $\kappa <1$, the refracted ray bends away from the normal. As a result, there is a critical value $\theta _c$ of the angle of incidence for which there is no value of $\theta _t$ unless $\sin \theta _i \leq \sin \theta _c = \kappa$. Thus for $\kappa < 1$ we impose the condition that

$$\begin{equation} x\cdot m \geq \kappa . {\tag{1.3}} \end{equation}$$

Similar physical restriction applies for $\kappa > 1$.

By using equation 1.1 or 1.2 we can define a map $T_{\mathcal{S}}: \Omega \to \Omega ^*$ that sends points in $\Omega$ (incident rays) to points in $\Omega ^*$ (reflected rays if we use 1.1 or refracted rays if we use 1.2) according to laws of geometric optics.

Generally, in an optical surface design problem the data consists of two bounded domains $\Omega$, $\Omega ^*$ contained in $\mathbb{R}^n$ and two probability measures $\rho _1, \rho _2$ defined on $\Omega$ and $\Omega ^*$ respectively, satisfying the energy conservation statement

$$\begin{equation} \rho _1(\Omega ) = \rho _2(\Omega ^*). \tag{1.4} \end{equation}$$

$\Omega$ and $\Omega ^*$ correspond to the directions of incident light and the directions of light after reflection or refraction respectively. Likewise $\rho _1$ and $\rho _2$ are the intensity of light from the source $\Omega$ and a desired illumination intensity on the target $\Omega ^*$, respectively. We regard $(\Omega , \rho _1)$ as a source pair and $(\Omega ^*, \rho _2)$ as a target pair.

The objective is then to find an optical surface (lens/mirror) $\mathcal{S}$ such that

$$\begin{equation} \rho _1(T_{\mathcal{S}}^{-1} (F)) =\rho _2(F) {\tag{1.5}} \end{equation}$$

for all $F \subset \Omega ^*$.

In practice, $\mathcal{S}$ is a reflecting mirror or refracting lens and for $m \in \Omega ^*$, $T_{\mathcal{S}}^{-1}(m)$ represents the set of all directions from the source that will, after reflection or refraction off of the lens $\mathcal{S}$, propagate in the direction of $m$. We use the notation $\mathcal{T}_{\mathcal{S}}$ for $T_{\mathcal{S}}^{-1}$. The map $\mathcal{T}_{\mathcal{S}}$ is called the ray-tracing map of $\mathcal{S}$.

The physical meaning of equation 1.5 is that the prescribed intensity of reflected or refracted light on the set $F \subset \Omega ^*$, which is $\rho _2(F)$, is equal to the intensity of light that comes from the source, which is $\rho _1(\mathcal{T}_{\mathcal{S}} (F))$. Depending on whether $\mathcal{S}$ is a reflector or a refractor, the design problem is called a reflector/refractor problem.

The models used to obtain numerical and analytical solutions for reflector/refractor problems are developed by using one of the following three approaches. One of the approaches is to use variational method where the problem is cast as an optimal mass transportation problem from $\rho _1$ to $\rho _2$ for an appropriate cost function and use Kantorovich duality, 20. Another method is to derive the second order nonlinear PDE of Monge–Ampère-type satisfied by the surface defining function from the energy conservation relation 1.5 between the light energy emitted by the source and the light energy received by the target and analyze the PDE, 11. A third approach is to exploit mainly the inherent geometric features of the problem and use a classical approach of Minkowski and Aleksandrov which was used in solving the Minkowski problem 16 of finding a closed convex surface whose Gaussian curvature is a given function of the exterior unit normal.

In the discussion below, we will focus on refractor problems and briefly describe the aforementioned approaches by using prototype refractor problems. We will state the problems explicitly and show how the available methods are used to solve the corresponding design problems. We will also mention problems that are worth studying in the future. We point out to the reader that unless mentioned otherwise, we consider $\kappa <1$.

2. Optimal Transport: The Far Field Refractor Problem

Several freeform lens design problems have been successfully modeled by using optimal transport framework on the sphere. See 41020. The models have also been adopted by the optics community. Among these problems is the far field refractor problem.

2.1. Statement of the far field refractor problem with point source

In this problem, we are given two domains $\Omega , \Omega ^*$ contained in the unit sphere $S^{n-1}$ of $\mathbb{R}^n$ (in practice $n=3$), $|\partial \Omega | =0$; a punctual source of light located at the origin $\mathcal{O}$, surrounded by media I, from which a monochromatic ray of light is issued in each direction $x \in \Omega$ with intensity density function given by $g(x)$ for $g \in L^1(\Omega ), g > 0$ a.e. on $\Omega$ and a prescribed intensity density distribution $f \in L^1(\Omega ^*)$ which satisfies

$$\begin{equation} \int _{\Omega } g(x) dx = \int _{\Omega ^*} f(x) dx . {\tag{2.1}} \end{equation}$$

The objective is to find a refracting surface parametrized as $\mathcal{S} = \{\rho (x)x : x \in \Omega \}, \rho >0$, between medium $I$ and medium $II$ such that all rays emitted from the point $\mathcal{O}$ with directions $x \in \Omega$ are transported via refraction by the surface into media $II$ with directions in $\Omega ^*$ in such a way that the prescribed illumination pattern $f \in L^1(\Omega ^*)$ is achieved on $\Omega ^*$.

In particular, if we let $d\rho _1 = g dx$ and $d\rho _2 = f dx$ where $dx$ is the surface measure on $S^{n-1}$, this means $\rho _1$ and $\rho _2$ satisfy the mass balance condition $\rho _1 (\Omega ) = \rho _2 (\Omega ^*)$ and a weak solution to the far field refractor problem is defined as an interface $\mathcal{S}$ so that the associated map $T_{\mathcal{S}}$ satisfies

$$\begin{equation} T_{\mathcal{S}}(\Omega ) \subset \Omega ^* \qquad \text{and} \qquad \rho _2 = (T_{\mathcal{S}})_\#\rho _1. {\tag{2.2}} \end{equation}$$

We recall for $F \subset \Omega ^*$ that the measure $(T_{\mathcal{S}})_\#\rho _1 (F) \coloneq \rho _1(T_{\mathcal{S}}^{-1}(F))$ is the push forward measure. Because of 1.3, we assume $x \cdot m \geq \kappa$ for $x \in \Omega$ and $m \in \Omega ^*$.

By using Snell’s law 1.2 it can be shown that the semi-ellipsoid of revolution (ellipse in 2D case) with one focus at $\mathcal{O}$, given in polar representation as (10)

$$\begin{equation*} E(m,b) = \left\{\dfrac{b}{1-\kappa m \cdot x}x: \ x \in S^{n-1},\ m \cdot x \geq \kappa \right\} \end{equation*}$$

for some $b>0$ and $m \in \Omega ^*$, has a uniform refracting property. This means that if $E(m,b)$ is an interface between media $I$ and $II$, all rays issued from $\mathcal{O}$ in a direction $x$ with $x\cdot m\geq \kappa$ will be refracted in the direction $m$ after hitting $E(m,b)$. Motivated by this, a solution for the far field refractor problem is sought among surfaces $\mathcal{S} = \{\rho (x)x : x \in \Omega \}$ which are supported from outside by semi-ellipsoids at each point. That is, for each point $\rho (x_o)x_o \in \mathcal{S}$ there exists a semi-ellipsoid $E(m,b)$ with $m \in \overline{\Omega ^*}$ such that $\rho (x) \leq \dfrac{b}{1-\kappa m \cdot x}$ with equality at $x_o$. We now discuss how optimal transport technique is used to find such a solution. First, a brief review of optimal transport.

2.2. Review of optimal transport problems

Let $X$ and $Y$ be compact subsets of $\mathbb{R}^n$ and $c: X \times Y \to [0, \infty ]$. We refer to $c$ as the cost function. Denote by $\mathcal{P}(X)$ and $\mathcal{P}(Y)$ the set of probability measures on $X$ and $Y$ respectively. Given $\mu \in \mathcal{P}(X)$ and $\nu \in \mathcal{P}(Y)$ a measurable map $T: X \to Y$ satisfying $T_{\#} \mu = \nu$ is called a transport map from $\mu$ to $\nu$. Monge’s formulation of the optimal transport problem is to minimize

$$\begin{equation} T \to \int _{X} c(x, T(x)) d\mu (x) {\tag{2.3}} \end{equation}$$

among all transport maps from $\mu$ to $\nu$. A minimizer of 2.3 is called an optimal transport map. A generalization due to Kantorovich 18 of Monge’s problem is stated as

$$\begin{equation} \inf _{\gamma \in \Gamma (\mu , \nu )} \int _{X \times Y} c(x, y) d \gamma (x, y) {\tag{2.4}} \end{equation}$$

among all transport plans $\gamma \in \Gamma (\mu , \nu )$ where

$$\begin{equation*} \Gamma (\mu , \nu ) : = \{\gamma \in \mathcal{P}(X \times Y) : (\pi _1)_{\#} \gamma = \mu , (\pi _2)_{\#} \gamma = \nu \}. \end{equation*}$$

Here $\pi _1: X \times Y \to X$ and $\pi _2: X \times Y \to Y$ are projection maps. If $c$ is lower semicontinuous and bounded from below, problem 2.4 has a minimizer in $\Gamma (\mu , \nu )$. A minimizer of 2.4 is called an optimal plan.

The $c$-transform of a given function $\phi : X \to [-\infty , \infty )$ is defined as $\phi ^c : Y \to \mathbb{R}$ where

$$\begin{equation*} \phi ^c(y) = \inf _{x \in X} c(x,y) - \phi (x). \end{equation*}$$

Also $\zeta ^{\bar{c}}: X \to \mathbb{R}$, the $\bar{c}$-transform of $\zeta :Y \to [-\infty , \infty )$, is defined by

$$\begin{equation*} \zeta ^{\bar{c}}(x) = \inf _{y \in Y} c(x,y) - \zeta (y). \end{equation*}$$

A function $\phi \in C(X)$ is called $c$-concave if $\phi = \zeta ^{\bar{c}}$ for some $\zeta$. We also define for any $x \in X$ the $c$-superdifferential $\phi$ at $x$ as

$$\begin{equation*} \partial _c \phi (x) = \{y \in Y: \phi (z) \leq c(z,y) - c(x,y) + \phi (x), \forall z \in X\}. \end{equation*}$$

Problem 2.4 is a linear optimization problem with convex constraints and thus, it admits a dual problem of maximizing

$$\begin{equation} \left\{ \int _X \phi (x) \, d\mu + \int _Y \psi (y) \, d\nu : \phi (x) + \psi (y) \leq c(x,y) \right\} {\tag{2.5}} \end{equation}$$

over $\phi \in L^1(X)$ and $\psi \in L^1(Y)$ 18.

One of the basic results in optimal transport theory is that if $c$ is continuous and $X$ and $Y$ are compact, then the dual problem 2.5 has a maximizer of the form $(\phi , \phi ^c)$ where the map $\phi$, which is referred to as Kantorovich potential, is a $c$-concave function. Furthermore, if $\partial _c \phi (x)$ is single valued for $\mu -a.e.\ x$, then $T: = \partial _c \phi$ induces the optimal plan $\gamma$ for 2.4. That is, $\gamma = (Id, T)_\# \mu$ minimizes 2.4. Consequently, $\partial _c \phi (x)$ is an optimal transport map from $\mu$ to $\nu$.

2.3. The far field refractor problem as optimal transport

To transform the far field refractor problem 2.2 to an optimal transport problem, we introduce the logarithmic cost function $c: \Omega \times \Omega ^* \to [0, \infty ]$ given by

$$\begin{equation} c(x,m) = -\log (1 - \kappa x \cdot m). {\tag{2.6}} \end{equation}$$

$c(x,m)$ is continuous and bounded from below since $x\cdot m \geq \kappa$. Thus, if we consider the dual problem 2.5 corresponding to $c(x,m)$ given as

$$\begin{equation} \sup _{u \in C(\bar{\Omega }), v \in C(\bar{\Omega }^*)} \left\{ \int _{\Omega } u(x) d\rho _1 + \int _{\Omega ^*} v(m) d\rho _2 : u(x) + v(m) \leq c(x,m)\right\}, {\tag{2.7}} \end{equation}$$

there exists 10 a $c$-concave $\phi$ such that $(\phi , \phi ^c) \in C(\bar{\Omega }) \times C(\bar{\Omega }^*)$ is a maximizer of 2.7. Let $\phi \in C(\Omega )$ be $c$-concave for the cost function $c(x,m) = -\log (1 - \kappa x \cdot m)$. If we set

$$\begin{equation} \rho (x) = e^{\phi (x)}, {\tag{2.8}} \end{equation}$$

then for each $x_o \in \Omega$ there exists $b>0$ and $m \in \Omega ^*$ such that

$$\begin{equation} e^{\phi (x)} \leq \dfrac{b}{1-\kappa m \cdot x} \,\, \forall x \in \Omega \,\, \text{ and } \,\, e^{\phi (x_o)} = \dfrac{b}{1-\kappa m \cdot x_o}. {\tag{2.9}} \end{equation}$$

That is, at each point $\rho (x_o) x_o$ the surface $\mathcal{S} = \{ \rho (x) x: x \in \Omega \}$ is supported from outside by a semi-ellipsoid $E(m,b)$ for some $b>0$ and $m \in \Omega ^*$.

If $\mathcal{S}$ is smooth, then $\mathcal{S}$ and $E(m,b)$ have the same normals at $\rho (x_o)x_o$. By the refraction property of $E(m,b)$ a ray emanating from $\mathcal{O}$ in a direction $x_o$ will be refracted off $\mathcal{S}$ in the direction $m$ by the surface $\mathcal{S} = \{x e^{\phi (x)}: x \in \Omega \}$. Therefore, $T_{\mathcal{S}}(x_o) = m$. Thus, $\mathcal{S}$ satisfies, $T_{\mathcal{S}}(\Omega ) \subset \Omega ^*$. Also, from 2.9 we get

$$\begin{equation*} \phi (x) \leq -\log (1 - \kappa m \cdot x) + \log (1 - \kappa m \cdot x_o) + \phi (x_o) \qquad \text{for all} \qquad x \in \Omega \end{equation*}$$

and hence

$$\begin{equation*} T_{\mathcal{S}}(x_o) = m \in \partial _c \phi (x_o). \end{equation*}$$

Furthermore if $c$ is given by 2.6, for any $c$-concave $\phi$ the $c$-superdifferential mapping $\partial _c \phi$ is single valued for $\rho _1-a. e.\ x$. Thus if $(\phi , \phi ^c)$ is a maximizer of 2.7, $\partial _c \phi$ is a measure preserving map from $\rho _1$ to $\rho _2$. Therefore $\partial _c \phi = T_{\mathcal{S}} \rho _1$-a.e., and more importantly $\rho _2 = (T_{\mathcal{S}})_\#\rho _1$ proving that $\mathcal{S} = \{x e^{\phi (x)}: x \in \Omega \}$ will be a solution of the far field refractor problem. From this we deduce that solving the far field refractor problem corresponds to finding a maximizer of the dual problem to the optimal transport problem 2.7. It is further shown in 10 that if $\mathcal{S} = \{x \rho (x): x \in \Omega \}$ and $\tilde{\mathcal{S}} = \{x \tilde{\rho }(x): x \in \Omega \}$ are two solutions of the far field refractor problem, then $\rho = C \tilde{\rho }$ for some constant C.

The optimal transport approach can also be used in the case of anisotropic materials where optical properties vary according to the direction of propagation of light 6. In this type of media modeling the design problem is more complicated due to issues such as birefringence, a situation where the incident rays may be refracted into two rays. It is also used to develop numerical methods to approximate solutions 4 to freeform design problems. Recently, a more efficient method called entropic regularization is applied to the the reflector problem 2 and it is likely that this approach can be applied to refractor problems.

Finally, it is worth noting that the far field problem has not been studied in its full generality using optimal transport method. In practice, if a ray of light hits an interface between two media which have different optical properties, the ray is partly reflected back and partly transmitted (Fig. 1). It is an open problem whether or not this energy imbalance between the source and target can be modeled as a variational problem by using the optimal transport approach.

3. Minkowski Approach: The Near Field Refractor Problem

The Minkowski problem involves finding a closed convex surface whose Gaussian curvature at a point with exterior normal $\xi$ is given by a continuous positive function $k(\xi )$ 16. To solve this, the approach was first to define a set function on $S^{n-1}$ by $\sigma (M) = \int _M \dfrac{d\omega (\xi )}{k(\xi )}$ where $d\omega$ is a surface element on the unit sphere, and consider the general problem of finding a convex hypersurface for which $\sigma$ is a surface area function. That is, to find a convex hypersurface $F$ such that for a given subset $M'$ of $F$, the area of $M'$ is equal to $\sigma (M)$ where $M' = \{x \in F: \text{the normal to $F$ at $x$ is in } M\}$. Then partition the sphere into smaller domains $g_k$ with area $\sigma _k$ and obtain an approximation $\sigma ' = \sum _{k=1}^\ell \sigma _k \delta _{\xi _k} ; \xi _k \in g_k$ to $\sigma$. Replacing $\sigma$ by $\sigma '$ gives a discrete problem of existence of convex polyhedron with faces of areas $\sigma _k$ and outward normal $\xi _k$. Let $P_\ell$ be a solution of the discrete problem. As $\ell \to \infty$ (or the diameter of $g_k$ to $0$), it is proved that the sequence $P_\ell$ of surfaces converges to a surface which solves the Minkowski problem. This classical approach of Minkowski is used to obtain iterative methods to solve not only various geometric optics problems related to both reflection 3 and refraction, but also to develop numerical methods to solve PDEs 15. The approach is also generalized to approximate solutions to semi-discrete optimal mass transport problems and generated Jacobian equations. See 1 and the reference therein. It is also one of the methods used by the optics community to design freeform optical surfaces 14. Below we exhibit how this approach is used in the near field refractor problem.

3.1. Statement of the near field refractor problem with point source

For this problem we are given $\Omega \subset S^{n-1}$, $\Omega ^* \subset \Sigma$, a hyperplane in $\mathbb{R}^n$ with $|\partial \Omega | =0$ and the origin $\mathcal{O} \notin \Omega ^*$. For each $x \in \Omega$, a light ray is issued in the direction of $x$ from a point source of light located at $\mathcal{O}$. The illumination intensity of the source is given by a density function $g \in L^1(\Omega )$, $g > 0$ a.e. A prescribed irradiance distribution on the target set $\Omega ^*$ is given by a nonnegative density function $f \in L^1(\Omega ^*)$. $f$ and $g$ satisfy the mass balance 2.1. We assume that $\Omega$ is surrounded by medium I and $\Omega ^*$ is surrounded by medium II. Additional constraints are imposed on $\Omega$ and $\Omega ^*$ to make the problem physically feasible 9.

The near field refractor problem involves finding an interface $\mathcal{S}$ between media $I$ and $II$, parametrized radially as $\mathcal{S} = \{\rho (x) x : x \in \overline{\Omega }\}$, that redirects each ray with direction $x\in \Omega$ by refraction into $\Omega ^*$ so that a prescribed irradiance distribution $f$ is obtained on $\Omega ^*$. More precisely we will require

$$\begin{equation} \int _{\mathcal{T}_{\mathcal{S}} (F)} g(x) \, d\sigma (x) = \int _F f(y) \, d\sigma (y) {\tag{3.1}} \end{equation}$$

holds for all Borel sets $F\subset \Omega ^*$.

Similar to the semi-ellipsoids in the far field problem, there are surfaces that have uniform refracting property in the near field case. The Descartes ovoid refracts all rays to a single point (Fig. 2). Let $P \in \Omega ^*$ and $\kappa |P| < b < |P|$. The Descartes ovoid with one focus at the origin $\mathcal{O}$ is given as

$$\begin{equation*} O_b(P) = \{h(x,P,b)x : x \in S^{n-1}, x\cdot P \geq b \} \end{equation*}$$

where

$$\begin{equation} h(x,P,b) = \dfrac{(b-\kappa ^2 x \cdot P)- \sqrt {(b-\kappa ^2 x \cdot P)^2 - (1-\kappa ^2) (b^2 - \kappa ^2 |P|^2)}}{1 - \kappa ^2}. {\tag{3.2}} \end{equation}$$

If the relative refractive index of the material inside the oval $\{h(x,P,b)x : x \in S^{n-1} \}$ to that of outside is $\kappa$, then by using 1.2 it can be shown that any ray emanating from the origin $\mathcal{O}$ and having direction $x\in S^{n-1}$ with $x\cdot P\geq b$ will pass through the point $P$ after refracting off of the ovoid $O_b(P)$.

Figure 2.

Refraction property of oval; $\kappa = 0.7$, $P = (10,0), b=8$.

The Descartes ovoids will be used as the building blocks of the near field refractors and we look for a solution for the near field refractor problem among surfaces which are supported from above by an ovoid at each point.

It is known that the near field refractor problem can’t be cast as an optimal mass transport problem. However, Minkowski’s approach can be used to solve the problem.

3.2. The near field refractor problem using Minkowski’s approach

Once again let $d\rho _1 = g d\sigma (x)$ and $d\rho _2 = f d\sigma (y)$. The road map is as follows. First partition $\Omega ^*$ and approximate the measure $\rho _2$ on the target by a sequence of discrete measures $\mu _{\ell }$ concentrated at finite points in $\Omega ^*$ and that converge to $\rho _2$ as $\ell \to \infty$. Then solve the problem when the target distribution is $\mu _{\ell }$. Finally let