1. INTRODUCTION

Lighting plays a crucial role in our current society, and since the introduction of light-emitting diodes (LEDs), the prevalence of beam-shaping optical elements has increased. This is partly because the sharp, point-like light from a bare LED package is typically considered undesirable and partly due to the increasing demand for aesthetic and personalized lighting, such as RGB LED lights. These optical elements are typically designed in an iterative and largely manual process, requiring significant experience and knowledge on the part of the optical designer, as well as considerable time ([1], Ch. 1.9). While a specular reflector can shape the light into a desired light distribution, it cannot necessarily reduce the sharpness of the light source since the mirrored surface may result in undesirable glare. Scattering elements, such as rough reflector surfaces exhibiting surface scattering or transmissive scattering elements, may be used to address this problem ([1], Ch. 1.8.4). Introducing scattering in the system generally means that some light control is lost, i.e., achieving the specified target cannot be guaranteed a priori. This work includes surface scattering in a consistent way in the existing framework for computing specular reflectors in the context of inverse freeform design to regain control over the light.

More precisely, the problem of directly computing an optical system given source and target distributions is often referred to as the inverse problem of illumination optics. Many methods of solving the inverse specular problem for reflectors and lenses have been developed over the last few decades, such as by solving a system of coupled ordinary differential equations (ODEs) in the case of rotationally or cylindrically symmetric systems [2]. For three-dimensional freeform optical surfaces, i.e., surfaces without any overall symmetry, a method that has proven successful is based on solving a Monge–Ampère equation [3–5].

While the specular inverse design problem is well researched, and a wide body of literature is concerned with predicting scattered light due to surface roughness within the physics community [6,7], and in the context of computer-generated imagery (GCI) [8], literature concerning the direct computation of scattering optical surfaces is scarce. The best reference we have found is Lin et al. [9], who designed a lens with a freeform scattering inner surface and a spherical outer one. Their approach represented the freeform surface by Bézier curves. The initial shape was iteratively modified to take into account the differences between the prescribed target distribution and the resulting raytraced distribution. While the aforementioned scattering models have proven very successful at their task of predicting scattered light, they are not always straightforward to apply in the context of computing freeform surfaces, as they often result in complicated expressions that are very difficult to invert. Additionally, many scattering models rely on wave optics, while freeform optical design via, e.g., a Monge–Ampère equation, usually considers the geometrical approximation of light where light propagates along straight lines called rays. Our goal with this work is to include scattering effects in freeform optical design and do so while remaining in the realms of geometrical optics. This approach allows us to utilize the vast amount of existing knowledge of freeform optical design.

As we showed in [10,11], the problem of computing two-dimensional, i.e., rotationally or cylindrically symmetric, reflectors with a scattering surface reduces to computing a deconvolution, followed by solving a specular reflector design problem. This paper will extend these results to compute three-dimensional freeform reflectors with scattering surfaces. To do so, we shall first find a mathematical relation between the light reflected from a perfectly smooth reflector with a specific shape and the scattered light from a reflector with the same shape made from a scattering material. We will show that this relation takes the form of a Fredholm integral equation of the first kind, and we will then show how we solved this integral relation to gain a suitable target function to use in the specular design problem. This approach is thus analogous to the two-dimensional one we presented in [10].

The paper is structured as follows. First, the proposed algorithm is summarized in Section 1.A. Then, the scattering model is derived in Section 2 based on concepts from optimal transport theory. Next, the freeform specular design problem is discussed in some detail in Section 3, followed by an outline of how we verified the scattering model in Section 4. Two numerical examples are shown in Section 5—the first showcases how we propose to use our model in a typical workflow, and the second shows how varying the amount of scattering influences the reflector shape. Finally, Section 6 concludes the paper with a summary and discussion of the topics presented herein, and some of our future ambitions.

A. Proposed Algorithm

This section briefly summarizes the proposed design algorithm. Given a desired far-field scattered light target distribution, a source distribution, and a probability density function (PDF) related to the scattering properties of the surface:

Assuming step 1 is well behaved, the above algorithm results in a reflector that transforms the source distribution into the original target distribution when taking scattering into account. The remainder of this paper is concerned with mathematically formulating each of these steps.

2. SCATTERING MODEL

This section treats the theoretical aspects needed to develop the scattering model and to apply it in the context of freeform reflector design.

A. Key Assumptions

We shall make several assumptions throughout our derivation of the Fredholm integral equation governing scattering in our model. The key assumptions are discussed here; additional assumptions will be introduced when they become relevant. The first assumption is that light scattering can be described using geometric optics by considering incoming, specularly reflected, and scattered light rays. Statistically, this is equivalent to the more physical notion of scattering whereby one incident direction yields multiple outgoing directions. Furthermore, light scattering is assumed to be fully elastic, i.e., the incident energy is scattered without absorption or other losses. The medium surrounding the reflector is also assumed to be lossless, and light is assumed to be scattered exclusively at the reflector surface. Finally, in this paper, we shall design only reflectors illuminated by zero-étendue parallel light with far-field targets. Note that the derivation of the scattering model is also valid for zero-étendue point sources since the scattering event occurs at the surface, irrespective of the system’s symmetry (or lack thereof).

B. Geometry

Suppose we have a Cartesian $xyz$-coordinate system in ${\mathbb{R}^3}$, with a parallel-ray source (henceforth parallel source) on a rectangular domain ${\cal S} = [{a_1},{a_2}] \times [{b_1},{b_2}] \subset {\mathbb{R}^2}$ in the $xy$ plane, centred around the origin ${\cal O}$; see Fig. 1. Rays emitted from the source (so-called source rays or incident rays) propagate in a fixed “upwards” direction, i.e., with a positive $z$ component, given by the unit vector ${{\hat {\textbf s}}}$ [hats (^) denote unit vectors throughout this paper]. For simplicity, we shall align the source rays with the $z$ axis, i.e., ${{\hat {\textbf s}}} \equiv {{\hat {\textbf{e}}}_z}{= (0,0,1)^ \intercal}$.

 

Fig. 1. Specular reflection from a smooth reflector.

Download Full Size | PDF

A specularly reflected ray (a so-called specular ray or reflected ray), i.e., one abiding by the familiar law of specular reflection, propagates along the unit vector

An off-specular ray leaving the surface (a so-called scattered ray) propagates along the unit vector

C. Model Derivation

We shall now consider the scattering model in detail. It is based on Monge’s formulation of the optimal transport problem in mathematics. The next few sections will show how this setup and subsequent analysis yield an expression for the scattered light in the form of a Fredholm integral equation of the first kind.

1. Mappings

Speaking in general terms, it can be shown that the optical map, i.e., the mapping that gives the specular direction $(\psi ,\chi)$ corresponding to a given incident direction $(\vartheta ,\varphi)$ parametrizing ${{\hat {\textbf s}}}(\vartheta ,\varphi)$, is injective for strictly convex mirrors (perfect specular reflectors) [12]. That is, the specular direction $(\psi ,\chi)$ is unique for any given incident direction. Let us denote this map by ${\textbf{m}}$, so that $(\psi ,\chi) = {\textbf{m}}(\vartheta ,\varphi)$.

Consider now a rough reflector, i.e., one where the resulting light is scattered into a direction $(\gamma ,\nu)$, typically different from $(\psi ,\chi)$. To relate the two directions, let us first define the elemental rotation matrices

Let $\alpha \in [0,\pi /2]$ and $\beta \in [0,2\pi]$ be the polar (incline) and azimuthal angles of the so-called cone vector

 

Fig. 2. Cone vector ${\hat {\textbf{c}}}$ traces a circle on the unit (hemi)sphere for fixed $\alpha$.

Download Full Size | PDF

Note that, by definition, ${\hat {\textbf{c}}}(\alpha ,\beta) = {{\textbf{R}}_z}(\beta){{\textbf{R}}_y}(\alpha){{\hat {\textbf{e}}}_z}$. Suppose we apply the rotation matrix ${{\textbf{R}}_y}(\psi)$ followed by ${{\textbf{R}}_z}(\chi)$ to ${\hat {\textbf{c}}}$. Then, for a fixed $\alpha \in [0,\pi /2]$, the resulting vector would trace a tilted cone coaxial with ${\hat {\textbf{t}}}$ by letting $0 \le \beta \le 2\pi$, with base radius $\sin (\alpha)$ and height $\cos (\alpha)$. This is what we want for our scattered vector ${\hat {\textbf{u}}}$, i.e.,

 

Fig. 3. Scattering by a rough reflector.

Download Full Size | PDF

Equating the representation of ${\hat {\textbf{u}}}$ in Eq. (7) with ${\hat {\textbf{u}}}(\gamma ,\nu)$ in Eq. (3) allows us to solve for $(\gamma ,\nu)$ for any given specular direction $(\psi ,\chi)$ and cone vector direction $(\alpha ,\beta)$. That is, we may find a so-called scattering map, ${\textbf{s}}$, such that ${\textbf{s}}(\psi ,\chi ,\alpha ,\beta) =\def\LDeqbreak{} (\gamma ,\nu)$. Suppose we instead want $(\alpha ,\beta)$ for some known $(\psi ,\chi)$ and $(\gamma ,\nu)$ pairs. This yields a third map, the so-called cone map, say ${\textbf{c}}$, such that ${\textbf{c}}(\psi ,\chi ,\gamma ,\nu) = (\alpha ,\beta)$. These relations are summarized in Fig. 4. We shall focus on the scattering part, i.e., finding ${\textbf{s}}(\psi ,\chi ,\alpha ,\beta)$ and ${\textbf{c}}(\psi ,\chi ,\gamma ,\nu)$.

 

Fig. 4. Relations between unit vectors and corresponding spherical coordinates.

Download Full Size | PDF

Finding the scattering map. We shall first find the scattering map, ${\textbf{s}}$, which returns the scattered direction $(\gamma ,\nu)$ of some specular direction $(\psi ,\chi)$ and cone direction $(\alpha ,\beta)$. Starting with $\gamma$, note that, by construction, $\cos (\gamma) = {u_3}$, where ${u_3}$ is the third component of ${\hat {\textbf{u}}}$. Thus, by evaluating Eq. (7), we find that

Finding the cone map. We shall now find the cone map, ${\textbf{c}}$, yielding $(\alpha ,\beta)$ for given directions $(\psi ,\chi)$ and $(\gamma ,\nu)$. Starting with $\alpha$, note that, by construction (recall Fig. 3), $\cos (\alpha) = {\hat {\textbf{t}}} \cdot {\hat {\textbf{u}}}$, so that

Finding $\beta$ requires significantly more effort. Theoretically, one could equate the two representations of ${\hat {\textbf{u}}}$ in Eqs. (3) and (7) and solve for $\beta$; in practice, however, this turns out to be very difficult. Since we shall assume rotational symmetry later in this paper, meaning the explicit expression for $\beta$ is no longer necessary, we only briefly summarize how it was derived below. We first considered a representation of ${\hat {\textbf{u}}}$ in terms of the stereographic components of the reflected vector ${\hat {\textbf{t}}}$ and the components of ${\hat {\textbf{c}}}$, i.e., ${\textbf{y}}({\hat {\textbf{t}}})$, ${c_1}$, ${c_2}$, and ${c_3}$, where ${\textbf{y}}$ is the 2-tuple associated with the unit vector ${\hat {\textbf{t}}}$ via stereographic projection from the north pole:

This allowed us to solve for ${c_1}$ and ${c_2}$ in terms of ${y_1}$ and ${y_2}$, or, via the stereographic projection, in terms of $\psi$ and $\chi$, as well as $\gamma$ and $\nu$. Then, we used the fact that, by construction, $\tan (\beta) = {c_2}/{c_1}$, and the scattering map in Eq. (11) to conclude that

2. Energy Balances

Let us introduce the light distributions associated with the source, specular and scattered light. The light source is parallel, meaning it is prescribed in the form of an exitance $[{\rm{W}} \cdot {{\rm{m}}^{- 2}}]$ denoted by $f(x,y) \gt 0$, $(x,y) \in {\cal S} \subset {\mathbb{R}^2}$, where ${\cal S} = [{a_1},{a_2}] \times [{b_1},{b_2}]$. Both the reflected light and scattered light may be described using intensity distributions [${\rm{W}} \cdot {\rm{s}}{{\rm{r}}^{- 1}}$] in the far field. We have:

Before discussing the freeform design problem, we shall formulate the relation between $g$ and $h$. To do so, let us first note that our assumptions regarding the absence of losses in the system lead to the following global energy balances:

Optimal transport. Suppose we fix $\psi = \Psi$ and $\chi = {\rm{X}}$ such that $(\Psi ,{\rm{X}}) \in {\cal T}$. Consider perfect specular reflection, i.e., reflection from a perfectly mirrored surface. In that case, $\alpha$ vanishes and $\beta$ is irrelevant, so that Eq. (7) gives ${\hat {\textbf{u}}} \equiv {\hat {\textbf{t}}}$. Then, $\gamma = \Gamma \equiv \Psi$ and $\nu = {\rm{N}} \equiv {\rm{X}}$, i.e., ${\textbf{s}}(\Psi ,{\rm{X}},\alpha ,\beta)$ is the identity map for all $\alpha \in [0,\pi /2]$ and $\beta \in [0,2\pi]$. This is true for any $\Psi$ and ${\rm{X}}$ so that we get the plots in Fig. 5.

 

Fig. 5. Specular maps $\psi \to \gamma$ and $\chi \to \nu$ with fixed points $\Psi \to \Gamma$ and ${\rm{X}} \to {\rm{N}}$.

Download Full Size | PDF

Suppose instead we have scattering from a rough surface; then $\alpha$ and $\beta$ are nonvanishing, and the simple one-to-one relationship schematically shown in Fig. 5 is replaced by a richer relationship. Schematically, we can imagine a broadening of the lines, indicating a probability to go in that direction—see Fig. 6. Here, $[{\Gamma _1},{\Gamma _2}] \times [{{\rm{N}}_1},{{\rm{N}}_2}] \subseteq {\cal U}$ indicates the nonzero region where the direction $(\Psi ,{\rm{X}}) \in {\cal T}$ is mapped. The relationship is more complex, but this is an intuitive starting point.

 

Fig. 6. Schematic scattering maps $\psi \to \gamma$ and $\chi \to \nu$ with fixed points $\Psi \to [{\Gamma _1},{\Gamma _2}]$ and ${\rm{X}} \to [{{\rm{N}}_1},{{\rm{N}}_2}]$.

Download Full Size | PDF

Let us now explore the connection between our approach of modeling scattering and optimal transport, particularly so-called Monge–Kantorovich problems ([13], Ch. 1). Let $\rho (\psi ,\chi ,\gamma ,\nu) \gt 0$, $(\psi ,\chi) \in {\cal T}$, $(\gamma ,\nu) \in {\cal U}$ represent a density with properties

Returning to the schematic scattering maps in Fig. 6, it seems reasonable to make the following ansatz. Let $p$ be a PDF on the unit sphere depicting the broadening of the lines. Then, the density $\rho$ that we shall pick is given by the product

 

Fig. 7. Scattering maps $\psi \to \gamma$ and $\chi \to \nu$ for various fixed values of $\chi$ and $\nu$ or $\psi$ and $\gamma$ from example #1 in Section 5.

Download Full Size | PDF

If we insert our choice of $\rho$ from Eq. (20) into the first relation of Eq. (18), meanwhile, and extend the limits to those of the unit sphere, we get, for all $\psi \in [0,\pi]$ and $\chi \in [0,2\pi]$,

Rotationally symmetric scattering. Suppose the PDF $p$ in Eq. (21) is rotationally symmetric, i.e., $p(\alpha ,\beta) = p(\alpha)$ for all $\beta \in [0,2\pi]$. In that case, Eq. (21) reduces to

Let us now return to the schematic scattering map in Fig. 6. In particular, compare the schematic versions to Fig. 7, which depicts the kernel $p$ from example #1 in Section 5. Since the kernel is a four-dimensional quantity via the mapping $\alpha (\psi ,\chi ,\gamma ,\nu)$, we consider only slices with two fixed angles—$\chi$ and $\nu$ or $\psi$ and $\gamma$. Focusing on the top row, it is clear that the mapping $\psi \to \chi$ does not change for slices where the azimuthal angles $\chi$ and $\nu$ are equal. However, the mapping significantly differs when $\chi \ne \nu$. In particular, values close to the poles are more likely to remain close to the poles and cannot readily reach, e.g., the equator. The situation is quite different for the mapping $\chi \to \nu$ (bottom row). First, notice that the map is naturally periodic at $\chi = 2\pi$ and $\nu = 2\pi$ since this is the period of the sphere. Next, note that these mappings do change for slices where the polar angles $\psi$ and $\gamma$ are equal. This is consistent with what is expected from the rotationally symmetric scattering PDF since we chose $p(\alpha)$ such that it is largest close to $\alpha = 0$ (i.e., $\psi = \gamma$) and then drops relatively rapidly to near-vanishing values—see example #1 in Section 5.

3. Unfolding the Fredholm Integral Equation

Suppose we want to solve the inverse problem; given a target intensity distribution $h$ and a scattering function $p$, can we compute the virtual specular distribution $g$? If the scattering equation had been a convolution integral, solving the inverse problem would be known as deconvolution. As they are Fredholm integrals, we shall refer to the process as unfolding for historical reasons—cf., e.g., [15]. Formally, there are constrained situations for which unique, closed-form, analytical solutions can be constructed when unfolding Fredholm integral equations ([16], Ch. 12), but we are interested in more general, numerical methods for obtaining an approximation of $g$.

We shall apply Richardson–Lucy deconvolution to the Fredholm integral problem. This method is based on maximum likelihood arguments; Richardson and Lucy independently developed a ratio deconvolution method due to a need for deblurring images from telescopes and in the context of fluorescence microscopy [17,18]. Specifically, Richardson used Bayesian statistics and assumed that a conditional probability caused the blurring, while Lucy considered maximizing the likelihood of the observed sample within the solution space. We shall not formally show that Richardson–Lucy deconvolution applies to unfolding our Fredholm integral equation. Still, the derivation by Lucy in [18] is so general that it is enough to assume that the blurring occurs via a Poisson process.

Before stating the final Richardson–Lucy expression, let us discretize the Fredholm integral equation in Eq. (21) [analogously for the rotationally symmetric case in Eq. (25)]. Fix a rectangular grid of ${N_1} \times {N_2}$ points, and let ${\textbf{h}}$ be the matrix representation of $h$ with components ${h_{\textit{ij}}} = h({\gamma _i},{\nu _j})\sin ({\gamma _i})$, where $i \in [1,{N_1}]$ and $j \in [1,{N_2}]$. Similarly, let ${\textbf{g}}$ be the matrix representation of $g$ such that ${g_{\textit{kl}}} = g({\psi _k},{\chi _l})\sin ({\psi _k})$, where $k \in [1,{N_1}]$ and $l \in [1,{N_2}]$. Finally, let ${\textbf{p}}$ be the tensor representation of $p$ such that $p_{\textit{ij}}^{\textit{kl}} = p({\textbf{c}}({\psi _k},{\chi _l},{\gamma _i},{\nu _j}))$, where $i,k \in [1,{N_1}]$ and $j,l \in [1,{N_2}]$. Then, the Fredholm integral can be written as

Let us now denote element-wise multiplication (i.e., Hadamard product) of two square matrices ${\textbf{A}}$ and ${\textbf{B}}$, fulfilling $\dim ({\textbf{A}}) = \dim ({\textbf{B}})$, as ${\textbf{A}} \odot {\textbf{B}}$. The resulting matrix has elements

It is worth noting that unfolding a Fredholm integral, much like deconvolution, is an ill-posed problem ([16], Ch. 12.12). In certain situations, one can show that deconvolution methods converge ([19], Ch. 5), but they are too restrictive to be of use to us and certainly do not apply to unfolding Fredholm integral equations using iterative deconvolution methods. We shall return to how well this approach works when discussing the numerical examples in Section 5.

3. FREEFORM SPECULAR REFLECTOR DESIGN

We shall now discuss how we computed the reflector surfaces that yield a desired target light distribution when light scattering is accounted for. First, Eq. (25) was discretized, then unfolded using the Richardson–Lucy method in Eq. (31), thus yielding a virtual specular target distribution ${{\textbf{g}}_{{\rm{uf}}}}$. Next, we computed a reflector that fulfills the resulting specular problem. To this end, we used the numerical Monge–Ampère solver first introduced in our group by Prins [3] and later expanded by Yadav [4] and Romijn [5]. The details are outside the scope of this paper, but a summary is given below.

A. Stereographic Coordinates

Because some equations become simpler to work within stereographic coordinates, these are used in the code for computing the reflector surfaces. In particular, we must transform our specular target intensity ${{\textbf{g}}_{{\rm{uf}}}}$ [${\rm{W}} \cdot {\rm{s}}{{\rm{r}}^{- 1}}$] into one defined in stereographic coordinates. Let us consider the continuous case with ${g_{{\rm{uf}}}}(\psi ,\chi)$. This section will omit the subscript $_{{\rm{uf}}}$ on $g$ to simplify the notation.

Let $f({\textbf{x}}) \gt 0$, ${\textbf{x}} \in {\cal S}$ be the source exitance distribution, and let $g(\psi ,\chi)$ be the specular target intensity (we shall return to the domain momentarily). Recall the geometry from Fig. 1, i.e., parallel rays leave the domain ${\cal S}$ in the $xy$ plane parallel to the positive $z$ axis and strike a reflector parametrized using some height function $z = u({\textbf{x}}) \gt 0$, ${\textbf{x}} \in {\cal S}$. Solving the inverse problem thus reduces to finding $u({\textbf{x}})$ such that $f$ is transformed into $g$.

Recall that ${\hat {\textbf{t}}}(\psi ,\chi)$ is the direction of the specular ray, and that ${\textbf{y}}$ is the 2-tuple stereographic representation of ${\hat {\textbf{t}}}$ defined in Eq. (13), where we have chosen stereographic projection from the north pole since the reflector surface is positioned above the parallel source, meaning the reflected rays typically travel “downwards,” i.e., in negative $z$ direction. Note that the stereographic projection from the north pole is undefined at the north pole itself, i.e., we consider ${t_3} \ne 1$ and $\psi \ne 0$. The corresponding inverse stereographic projection is

B. Energy Conservation

Let $\tilde g$ be the stereographic representation of $g(\psi ,\chi)$ such that $\tilde g({\textbf{y}}) = g(\psi ({\textbf{y}}),\chi ({\textbf{y}}))$. Let ${\cal A} \subseteq {\cal S}$ be a (sub)set of the source domain. Local energy conservation in the far-field approximation then states

Let the optical map in stereographic coordinates be ${\textbf{y}} = \widetilde {\textbf{m}}({\textbf{x}})$ [recall that $(\psi ,\chi) = {\textbf{m}}({\textbf{x}})$ is the optical map in spherical coordinates]. Then, Eq. (34) becomes (after substitution and transformation to integration over ${\textbf{x}}$)

C. Numerical Solution to Monge–Ampère Equation

Solving the standard Monge–Ampère equation in Eq. (38) is nontrivial, and thus a numerical least-squares approach was chosen. We need to venture further outside the scope of this paper to describe the method in detail. Thus, we point the reader to the works of Prins, Yadav, and Romijn [3–5]. Note, however, that we shall always compute the strictly convex solution, such that ${\textbf{m}}$ and $\widetilde {\textbf{m}}$ are injective mappings.

4. VERIFICATION

This section shows how we numerically verified our scattering model in Section 5. In particular, once the reflectors had been computed in the manner described in Section 3, they were raytraced, and the resulting distributions were then compared to the ones predicted by our model.

A. Raytracing

To this end, we wrote a custom raytracer that directly implements the model of scattering presented in Section 2. This approach was chosen instead of using pre-existing raytracing software such as LightTools to completely control the scattering behavior so that the model could be reliably verified. We note that our “raytracer” is much simpler than any commercial tool since we need it to work only for a parallel source. This allows us to not trace rays in the traditional sense—we never compute ray–triangle intersections, for instance—but rather determine the direction of the rays using the relations for ${\hat {\textbf{t}}}$ and ${\hat {\textbf{u}}}$. Since our targets are in the far field, it is enough to compute the spherical angles of ${\hat {\textbf{t}}}$ and ${\hat {\textbf{u}}}$ to compare the reflector’s performance versus what we prescribed.

1. Implementation

The raytracer was implemented in MATLAB, and it works as follows. First, the normals of the reflector are computed for each grid point on the rectangular grid using MATLAB’s surfnorm routine. Next, an arbitrary point $(x,y)$ in the source domain ${\cal S}$ is sampled using MATLAB’s rand command (for simplicity, we always use a constant source). The normal ${\hat {\textbf{n}}}(x,y)$ is then found using MATLAB’s interp2 routine with piecewise linear interpolation, i.e., each component of the normal vector is assumed to change linearly between the closest known normals on the initial rectangular grid. The reflected direction ${\hat {\textbf{t}}}(\psi ,\chi)$ associated with this source ray is computed using the vectorial LoR, Eq. (2)—note that ${{\hat {\textbf s}}} \equiv {{\hat {\textbf{e}}}_z}$ for the parallel source. Next, the scattered ray is computed by applying Eq. (7) using $\psi$ and $\chi$ from ${\hat {\textbf{t}}}$, and $\alpha$ and $\beta$ by sampling from the appropriate PDF $p$—see the next paragraph for details on this sampling.

We discretized the domains ${\cal S}$, ${\cal T}$, and ${\cal U}$ to collect the source, specular, and scattered rays, forming so-called “bins.” To determine the bin corresponding to a given direction, we use MATLAB’s dsearchn nearest-point search routine, and then the number of rays in the returned bin is incremented.

Once the desired number of rays has been traced, the ray count per bin is converted into an exitance or intensity, such that it may be compared to $f$, $g$, or $h$. For instance, suppose we have an ${N_1} \times {N_2}$ grid of bins. Then, the exitance of the source is estimated using ($i \in [1,{N_1}]$, $j \in [1,{N_2}]$)

Sampling of $\alpha$ and $\beta$. We shall now consider the problem of sampling $\alpha$ and $\beta$ in our raytracer. We have already mentioned that $\beta$ will be sampled uniformly on $[0,2\pi]$, so we need to consider only how $\alpha$ is sampled. For instance, suppose we would like a “rotationally symmetric Gaussian on the sphere.” This is a vague definition that can be interpreted in many ways. For example, we could pick $\alpha$ from a regular one-dimensional Gaussian via a process known as inverse transform sampling and then pick $\beta$ uniformly, or we could use a generalized distribution for picking normally distributed points on a sphere, such as the Kent distribution used in geology and bioinformatics [21].

We have instead chosen the following approach. First, we pick two independent normally distributed variables ${q_1} \sim {\cal N}(0,\sigma)$ and ${q_2} \sim {\cal N}(0,\sigma)$, where ${\cal N}(\mu ,\sigma)$ is the normal distribution with mean $\mu$ and standard deviation $\sigma$. The point ${\textbf{q}}: = ({q_1},{q_2}{)^ \intercal} \in {\mathbb{R}^2}$ in the $xy$ plane is picked from a rotationally symmetric two-dimensional normal distribution. Applying inverse stereographic projection from the south pole to ${\textbf{q}}$ thus yields a point on the unit sphere, representing the direction of the cone vector ${\hat {\textbf{c}}}$. That is,

B. RMS Error

The root mean square (RMS) error will quantify the error between the raytraced and exact distributions. For a discretization grid of ${N_1} \times {N_2}$, i.e., ${N_1}$ polar angles and ${N_2}$ azimuthal angles, the RMS error between ${\textbf{h}}$ and the raytraced ${{\textbf{h}}^*}$ is defined as

5. NUMERICAL EXAMPLES

This section discusses two numerical examples to showcase the design procedure outlined in this paper and the effect that varying amounts of surface scattering have on the shape of the computed freeform reflectors. Let

A. Example #1: Overlapping Gaussians

The first example we considered is outlined below, where $p(\alpha ;\sigma)$ can be found in Eq. (43).

The prescribed distributions are shown in Fig. 8, where we opted to plot $p$ on the unit sphere to facilitate comparisons to Fig. 2. Note that the most likely locations for the cone vector are close to the $z$ axis, i.e., relatively small-angle scattering, and we can be confident that $\alpha \lt \pi /2$. In addition to the false-color plots, we have also “sliced” each distribution along the red lines with corresponding plots to the right of the accompanying false-color plots. Note that we choose $f$ and $h$ such that (energy conservation)

 

Fig. 8. Prescribed distributions in example #1; ${128^2}$ sample points.

Download Full Size | PDF

Next, we computed the “unfolded” distribution ${g_{{\rm{uf}}}}$ using 1000 Richardson–Lucy iterations, i.e., by applying Eq. (31) 1000 times. This yielded the distributions in Fig. 9. Note that we have nothing to directly compare ${g_{{\rm{uf}}}}$ to since the exact solution is unknown for this problem. Thus, we computed the so-called “refolded” ${{\textbf{h}}_{{\rm{rf}}}}: = {\textbf{p}}{{\textbf{g}}_{{\rm{uf}}}}$, representing the scattered distribution that would occur from a reflector designed using ${g_{{\rm{uf}}}}$. This is shown in Fig. 10. As we can see, $h$ and ${h_{{\rm{rf}}}}$ are very similar, meaning ${g_{{\rm{uf}}}}$ is a good representation of the “true” $g$—at least in the sense that the predicted scattered distribution from the reflector will be close to the prescribed target distribution. This indicates that the Richardson–Lucy deconvolution method works well for the more general problem of unfolding our Fredholm integral.

 

Fig. 9. Unfolded specular target distribution ${g_{{\rm{uf}}}}$ in example #1; ${128^2}$ sample points, 1000 Richardson–Lucy iterations.

Download Full Size | PDF

 

Fig. 10. Refolded scattered target distribution ${h_{{\rm{rf}}}}$ in example #1; ${128^2}$ sample points.

Download Full Size | PDF

Next, we computed the specular reflector that achieves ${g_{{\rm{uf}}}}$ in the far field, given $f$ as the source. Because of the way the least-squares algorithm works, we must specify a boundary in the target domain and not have values too close to zero of ${\tilde g_{{\rm{uf}}}}$ within this boundary—recall Eq. (38). Since our initial scattered target distribution consisted of overlapping Gaussians, the support of ${g_{{\rm{uf}}}}$ is not finite. We thus chose some $\epsilon$ as a cutoff and found the boundary outlining the specular target domain by

 

Fig. 11. Specular target together with the support of ${g_{{\rm{uf}}}}$ in example #1; ${128^2}$ sample points.

Download Full Size | PDF

Since we altered ${g_{{\rm{uf}}}}$ by giving it a finite domain and renormalizing, we must naturally update our predicted scattered light distribution accordingly, i.e., apply the Fredholm integral equation again, with the new ${g_{{\rm{uf}}}}$ as the specular target distribution. The result is shown in Fig. 12, where we can see that the final predicted ${h_{{\rm{rf}}}}$ is very similar to the original $h$, but with slightly increased values close to the center and somewhat decreased values further out due to the cut tails of the specular target ${g_{{\rm{uf}}}}$.

 

Fig. 12. Final predicted scattered distribution in example #1; ${128^2}$ sample points.

Download Full Size | PDF

Next, the reflector was computed using the least-squares solver, and then the custom raytracer was employed to verify the shape of the surface and our scattering model. The results of the ray trace after ${10^7}$ rays are shown in Fig. 13, where we see that the source sampling is homogeneous and that the scattered light distribution is very close to our prediction, as is also confirmed by the RMS error essentially following the expected $N_r^{- 1/2}$ trend, where ${N_r}$ is the number of rays traced ([20], p. 9). The specular distribution does deviate from the prescribed target in some places, which is best demonstrated by the slices—see especially the data points close to the boundary of ${\cal T}$ and those close to the peaks. These discrepancies come from the numerical least-squares solver used to compute the reflector surface. They are presumably the source of the slight deviation from the theoretical $N_r^{- 1/2}$ slope of the RMS error, too. The results could perhaps be improved by using a finer grid. Crucially, however, the discrepancies do not significantly affect the scattered light, which is the main topic of concern here—it is still clear that our predictions for the scattered light align very well with what is obtained from the raytracer.

 

Fig. 13. Raytraced distributions and the reflector in example #1; ${128^2}$ sample points, ${10^7}$ rays traced.

Download Full Size | PDF

B. Example #2: Varying Amounts of Scattering

In our second example, we wanted to visualize and quantify the differences in the reflector shape due to varying the amount of surface scattering. The problem is outlined below, where we again enforced energy conservation such that the integrals over the source and target distributions were unity. The surface scattering function $p$ can be found in Eq. (43), with varying $\sigma$, i.e., varying amounts of scattering in the system. Here, $\sigma = 0$ signifies a smooth, specular reflector.

The PDFs $p$ with nonzero values of $\sigma$ and the associated specular target distributions are shown in Fig. 14, together with the boundary of ${\cal T}$, found by fixing $\epsilon = 0.1\;\max ({g_{{\rm{uf}}}})$. Successively increasing $\sigma$ shows a “sharpening” of the target; the target domain shrinks, and the maximum value increases.

 

Fig. 14. Probability density functions and specular target distributions in example #2; $\sigma \in \{0.025,0.05,0.075,0.1\}$ from left to right; ${64^2}$ sample points.

Download Full Size | PDF

Figure 15 shows the effect that scattering has on the reflector, where the base reflector was taken as the specular reflector achieving $h$ given $f$, i.e., with $\sigma = 0$ so that $p$ is a delta function and $g \equiv h$. Successive reflectors have increasing values of $\sigma$, associated with more and more scattering up to $\sigma = 0.1$. As expected, more scattering requires more modification of the reflector versus the base one, and we see variations in height up to a few percent of the size of the reflectors, which is consistent with previous observations we made in the two-dimensional case [11]. As noted there, variations of this order of magnitude are typically considered manufacturable.

 

Fig. 15. Slices along the indicated lines for reflectors with varying amounts of surface scattering, all fulfilling the problem in example #2; ${32^2}$ sample points.

Download Full Size | PDF

6. CONCLUSIONS AND DISCUSSION

We have developed a novel approach to computing freeform reflectors with scattering surfaces inspired by optimal transport. Our method involves using a density function with specific properties, which results in a Fredholm integral equation of the first kind. This equation provides information about scattered light in the far field based on a PDF that defines the surface’s scattering properties and the specular distribution. We can create a virtual specular target distribution by unfolding (i.e., solving) the Fredholm integral, in a process analogous to deconvolution. This virtual distribution can then be used as a target when solving the inverse problem of illumination optics, i.e., computing the reflector. This approach ensures that the prescribed target is achieved when considering surface scattering.

As a result, the process of designing freeform reflectors with scattering surfaces becomes a two-step process. We first modify the target distribution by unfolding the Fredholm integral equation that governs scattering in our model. Then, we use existing methods to compute the associated specular reflector. The addition of scattering enables new design possibilities for optical engineers, and the relative simplicity of our approach means that it can be readily implemented in a real workflow. As always, simplicity comes with a price—for our model to be valid as formulated herein, the scattering PDF must hold everywhere on the surface, i.e., local changes in scattering properties are not considered. In contrast to the approach of Lin et al. [9], we focused on computing freeform scattering reflectors as opposed to freeform lenses. However, this is not an inherent restriction in our model, meaning extending it to lenses exhibiting surface scattering should be possible.

Our future goal is to expand on the approach we introduced in [11] by applying it to three dimensions. In that work, we utilized microfacets to model the rough surface that causes light scattering. Furthermore, one way to enhance the applicability of our model is by removing the current limitation of isotropic surfaces.