1. INTRODUCTION

The inverse problem of light transport within the field of computational illumination optics is concerned with computing an optical system—typically lenses or reflectors—such that a given source light distribution is transformed into a desired target distribution [1–5]. A common situation involves shaping light from a point source since this is the approximate distribution of a light-emitting diode (LED). Romijn et al. extensively studied shaping point source illumination using reflectors and lenses by constructing numerical solutions to the so-called generalized Monge–Ampère equations [4,6]. For cylindrically and rotationally symmetric reflector problems, we can restrict our attention to two dimensions and solve two coupled ordinary differential equations (ODEs) to compute the specular surface [7,8].

Since the introduction of LED light sources, there is now significant demand for aesthetic lighting in terms of color and light shaping. Diffusive media or scattering surfaces can be utilized to homogenize light from an LED source or hide the LED itself ([9], p. 24). Additionally, no actual reflector surface is a perfect mirror—they all exhibit some form of light scattering due to surface defects introduced in the manufacturing process ([9,10], p. 24). Since the aforementioned methods of computing reflectors—i.e., Monge–Ampère equations and coupled ODEs—only consider specular reflection from a perfect mirror, the effect of light scattering due to surface roughness is absent. As such, real reflectors manufactured according to these designs may result in a broadened outgoing light distribution and blurring of finer features due to surface defects ([10], Sec. 1.8.4), ([11], Ch. 10). The present work aims to integrate surface light scattering into the design process of two-dimensional freeform reflectors in a mathematical way, thus allowing scattering surfaces to be utilized in their designs. Literature concerning the unification of scattering and freeform optics is minimal. The best reference we are aware of is Lin et al., who designed a lens with a freeform scattering inner surface and a spherical outer one [12]. The freeform surface was initialized using a Bézier curves representation and subsequently iteratively modified to take into account the differences between the prescribed target distribution and the resulting raytraced distribution.

In contrast, we will focus on computing two-dimensional freeform reflectors with a scattering surface using a direct inverse method. Specifically, we will derive a surface light scattering model by endowing a smooth curve (representing the surface) with microscopic roughness using so-called microfacets and consider how the outgoing light distribution changes by their orientations. The resulting expression for the scattered light is a convolution integral between a probability density function (PDF) dictating the microfacet orientations (and is, hence, related to the surface roughness) and a so-called virtual specular light distribution representing the reflected light from a perfectly smooth reflector. Thus, assuming the PDF can be estimated a priori, deconvolution can be utilized to compute a virtual specular light distribution, which can subsequently be inserted as the target distribution in the inverse two-dimensional specular reflector problem solved by Maes [8].

The structure of the manuscript is as follows. First, the theory needed to develop the surface scattering model is introduced (zero-ètendue sources—Section 2.C.1, the bidirectional reflectance distribution function (BRDF)—Section 2.C.2, and the microfacets—Section 2.C.3). Then, the surface scattering model is derived for parallel-ray and point sources in two dimensions (Section 2.D.1), and the necessary two-dimensional specular design ODEs are derived (Sections 2.D.2 and 2.D.3). Next, several aspects of the model are verified using raytracing by considering three numerical examples (Section 3), and finally some conclusions and proposed extensions are presented in Section 4.

2. SCATTERING MODEL AND FREEFORM REFLECTOR DESIGN

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

A. Key Assumptions

We have made several assumptions throughout our derivation of the convolution integral governing scattering in our model. The key assumptions are discussed here; any 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 rays. Statistically, this is equivalent to the more physical notion of one incident direction yielding multiple outgoing directions. Additionally, light scattering is assumed to be fully elastic, i.e., the incident energy is scattered without absorption or other losses. The medium surrounding our reflector is also assumed to be lossless, and light is assumed to be scattered exclusively due to surface roughness. Finally, the addition of scattering in the system is assumed not to break the two-dimensional nature of the problem, i.e., we consider in-plane scattering.

We postulate that the assumption of in-plane scattering holds in situations where the overall two-dimensional symmetry remains, such as illumination by a line source onto an extruded reflector. The fact that the specular reflector problem for a cylindrically symmetric system can be considered in cross section has been shown in, e.g., [7,8]. Note that the line source is parallel to the extrusion direction and suppose that the process of extrusion results in grooves in the reflector surface. These grooves are then also parallel to the extrusion direction, so the symmetry is preserved, and we postulate that we may study a two-dimensional cross section of the three-dimensional problem.

B. Geometry

The two-dimensional geometry of scattering used throughout this manuscript is outlined in Fig. 1.

 

Fig. 1. Geometry of scattering; ${\vartheta _r} = - {\vartheta _i},{\vartheta _o} \lt 0$, and ${\vartheta _i} \gt 0$.

Download Full Size | PDF

Here, we see an incident ray along direction ${\hat{\textbf s}}$ striking the surface at a point ${\cal P}$, resulting in a specular ray along ${\hat{\textbf t}}$ for perfect specular reflection or a scattered ray along ${\hat{\textbf u}}$ for rough surfaces. The directions ${\hat{\textbf s}}$, ${\hat{\textbf t}}$, and ${\hat{\textbf u}}$ are given by angles ${\vartheta _i} \in (- \pi /2,\pi /2)$, ${\vartheta _r} = - {\vartheta _i}$, and ${\vartheta _o} \in (- \pi /2,\pi /2)$, respectively, i.e.,

Note that hats $(\;\hat {}\;)$ denote unit vectors throughout this manuscript.

C. Scattering Model

Having introduced the geometry, let us now consider how the model is constructed. We will first introduce the concept of zero-étendue sources, followed by the BRDF, and finally the microfacets used to model surface roughness.

1. Zero-Étendue Sources

The concept of étendue is easiest to consider in phase space—the space in which all possible states of the system are represented. The phase space of our two-dimensional geometry consists of two dimensions—one spatial coordinate $x$ representing the position and one coordinate representing the momentum $p: = n\sin (\vartheta)$ of the ray at $x$ in direction $\vartheta$, where $n$ is the refractive index. We will always take $n$ as unity so that $p \equiv \sin (\vartheta)$.

The étendue in our two-dimensional phase space is given by {cf. [13], Eq. (5.50), p. 89}

 

Fig. 2. Illustration of a parallel source of length ${\ell _s}$ emitting rays in direction $\vartheta = {\vartheta _s}$ (left) and a point source at $x = {x_s}$ emitting rays in directions $\vartheta \in [{\vartheta _1},{\vartheta _2}]$ (right); ${\vartheta _1}$, ${\vartheta _s} \lt 0$ and ${\vartheta _2} \gt 0$.

Download Full Size | PDF

 

Fig. 3. Nonvanishing regions of phase space for the parallel (left) and point (right) sources in Fig. 2.

Download Full Size | PDF

Definition 1. The radiance of a two-dimensional parallel source emitting rays in direction ${\vartheta _s}$ is given by

Definition 2. The radiance of a two-dimensional point source at position ${x_s}$ is given by

The nonvanishing regions of phase space for the parallel and point sources in Fig. 2 are shown in Fig. 3, i.e., a horizontal and vertical line segment, respectively. Returning to the definition of étendue, Eq. (2), it is clear that it vanishes for both these sources since they form line segments of zero width in phase space—the parallel source has vanishing angular extent, and the point source has vanishing spatial extent. Hence, they are so-called zero-étendue sources—the only type of light sources we will consider in this manuscript.

2. Bidirectional Reflectance Distribution Function

The BRDF $B({\vartheta _i},{\vartheta _o})$ [${{\rm rad}^{- 1}}$] was first formulated by Nicodemus in 1965 [14], and it describes the way light from direction ${\vartheta _i}$ is scattered from the surface into direction ${\vartheta _o}$ [15] (recall Fig. 1). To see how, let ${L_i}({\vartheta _i})$ [${\rm W} \cdot {{\rm rad}^{- 1}} \cdot {{\rm m}^{- 1}}$] be the incoming radiance from direction ${\vartheta _i}$, and let ${L_o}({\vartheta _o})$ [${\rm W} \cdot {{\rm rad}^{- 1}} \cdot {{\rm m}^{- 1}}$] be the outgoing (i.e., scattered) radiance in direction ${\vartheta _o}$. Then, the so-called scattering equation, yielding the outgoing radiance in direction ${\vartheta _o}$, may be formulated as {cf. [14], Eq. (7)}

It can be shown that any physical BRDF must satisfy the following properties [16]:

Specular reflection as a BRDF. For example, let us consider the BRDF of specular reflection from a (locally) flat surface—see Fig. 4. The familiar law of reflection yields ${\vartheta _o} \equiv {\vartheta _r} = - {\vartheta _i}$. This leads to the following definition for the BRDF of specular reflection.

 

Fig. 4. Perfect specular reflection from a flat surface; ${\vartheta _o} = - {\vartheta _i} \lt 0$ and ${\vartheta _i} \gt 0$.

Download Full Size | PDF

Definition 3. The BRDF of specular reflection is given by

This BRDF fulfils:

Additionally, the scattering equation, Eq. (5), yields

Note that we are keeping the cosine terms in this relation for later, but we could write ${L_o}({\vartheta _o}) = {L_i}({\vartheta _i})$ here. That is, the radiance in the outgoing direction ${\vartheta _o}$ is exactly the radiance from the incoming direction ${\vartheta _i}$, which is what we expect for specular reflection.

3. Microfacets

Having treated specular reflection, we are now ready to develop a simple surface roughness model based on the concept of microfacets, i.e., small, tilted, specular sections superimposed on the macroscopic reflector. By describing the orientations of the microfacets along the reflector surface using some suitable PDF, a realistic model for surface roughness may be constructed, as has been demonstrated by its extensive use in computer-generated imagery (CGI) ([17], Ch. 8.4). Of particular interest to us is the microfacet model developed by Torrance and Sparrow in 1967 [18] and later applied to CGI by Cook and Torrance in 1982 [19]. We have opted to develop a simpler model inspired by the one introduced by Torrance and Sparrow (in particular, we did not consider effects such as shadowing and masking) since our primary goal is to unify surface scattering with inverse design, starting with as few complexities as possible while remaining realistic. Note that the simplicity of our surface roughness model means that it is not straightforward to directly relate it to, e.g., the root mean square (RMS) roughness of a real surface. In CGI, a commonly used microfacet distribution function based on the assumption of a Gaussian distribution of microfacet slopes is the Beckmann–Spizzichino model ([20], Ch. 5). In contrast to our simple model, this model contains a term proportional to the RMS slope of the microfacets {[17], Eq. (8.9)}.

To describe the microfacets, consider zooming in sufficiently close so we can observe a locally flat section of the rough reflector. Next, rotate the view such that the unit normal of the reflector, ${\hat{\textbf n}}$, aligns with the $y$ axis; this yields the view in Fig. 4 and defines our coordinate system with its origin at ${\cal P}$. Finally, consider tilting this locally flat reflector piece by some $\eta \in (- \pi /2,\pi /2)$, while keeping the angle of the incident light constant with respect to the original unit normal ${\hat{\textbf n}}$—see Fig. 5, where $\vartheta _i^m$ and $\vartheta _o^m$ are the incident and outgoing directions with respect to the microfacet normal ${{\hat{\textbf n}}^m}$, respectively. Mathematically, we define a rotation matrix as

 

Fig. 5. Geometry of a microfacet (solid black line) versus the macroscopic surface (dotted black line); ${\vartheta _o},\vartheta _o^m \lt 0$ and $\eta ,{\vartheta _i},\vartheta _i^m \gt 0$.

Download Full Size | PDF

Since the microfacet is specular, the law of reflection yields $\vartheta _o^m \equiv \vartheta _r^m = - \vartheta _i^m$. Using this, together with the geometry in Fig. 5, we get

That is, the outgoing direction with respect to the original macroscopic normal is the negative incident direction, plus 2 times the perturbation of the normal. Inserting this expression for ${\vartheta _i}$ in Eq. (12) yields

This relation gives the outgoing radiance ${L_o}$ in terms of the incoming radiance ${L_i}$ for a microfacet with orientation $\eta$ with respect to the macroscopic surface normal.

Distributed orientations. To use microfacets to describe scattering from a surface, rather than fixing $\eta$, we will consider that it is sampled from some distribution related to the roughness properties of the surface. That is, we do not know a priori which microfacet orientation we strike, but rather we wish to describe the distribution of possible orientations using some PDF $\tilde p(\eta)$ with support ${\rm supp}(\tilde p) = [{\eta _1},{\eta _2}] \subseteq (- \pi /2,\pi /2)$, i.e., $\int_{{\eta _1}}^{{\eta _2}} \tilde p(\eta) {\rm d}\eta = 1$. Note that, while the support of $\tilde p$ can formally extend to $(- \pi /2,\pi /2)$, we will typically consider situations where the significant contributions are concentrated in a much more narrow region centered around $\eta = 0$. This shape is consistent with measured BRDFs from machined surfaces ([11], Ch. 4).

To get the outgoing radiance, we now multiply the right-hand side of Eq. (15) by $\tilde p(\eta)$ and integrate over all valid microfacet orientations, leading to the below definition.

Definition 4. Let $\eta \in [{\eta _1},{\eta _2}] \subseteq (- \pi /2,\pi /2)$ be the microfacet orientation with respect to the macroscopic normal ${\hat{\textbf n}}$, sampled from the PDF $\tilde p(\eta)$. Let ${L_o}$ be the outgoing (i.e., scattered) radiance in direction ${\vartheta _o}$ due to an incoming radiance ${L_i}$ from direction ${\vartheta _i}$. Then,

By using $\eta \equiv ({\vartheta _i} + {\vartheta _o})/2$—recall Eq. (14)—and changing the integration variable to ${\vartheta _i}$, we get

Since ${\vartheta _o} \in (- \pi /2,\pi /2)$, and ${\rm supp}(\tilde p) = [{\eta _1},{\eta _2}] \subseteq (- \pi /2,\pi /2)$, we can rewrite this integral as

Note that the above expression does not consider whether a ray can physically reach ${\vartheta _o}$ from ${\vartheta _i}$. As an extreme example, consider near-grazing incident light, where any nonzero perturbation $\eta$ of the normal could mean the incident ray strikes the back of the microfacet. Implicit in the derivation of this model is that such unphysical situations do not occur.

Next, let

Before proceeding, let us briefly verify that $p$ is also a PDF. Let $\alpha : = {\vartheta _i} + {\vartheta _o} = 2\eta$, then

D. Freeform Reflector Design

Thus far, we have derived the scattered radiance [${\rm W} \cdot {{\rm rad}^{- 1}} \cdot {{\rm m}^{- 1}}$] from a single point—Eq. (21). Let us now consider what occurs when the point is located on a curve that traces a reflector. We will show how this situation yields a convolution integral for the outgoing intensity [${\rm W} \cdot {{\rm rad}^{- 1}}$], which can be used to design reflectors via an inverse specular problem by using deconvolution. We will then derive the relations needed to compute the reflectors for both parallel and point sources. More details about two-dimensional specular freeform reflector design can be found in Maes’s thesis ([21], Chs. 1–6).

1. Intensities

Let $s \in [0,\ell]$, parametrizing a curve of length $\ell$ be the arc length. The point-wise microfacet scattering equation, Eq. (21), then becomes

For a perfectly smooth reflector, an incident ray is specularly reflected from direction ${\vartheta _i}$ into direction ${\vartheta _r} = - {\vartheta _i}$ at every point $s$. Defining a virtual specularly reflected radiance ${L_r}(s,{\vartheta _r}): = {L_i}(s, - {\vartheta _r})$, allows us to write the above relation as

This convolution integral allows the interpretation that an incident ray is first reflected into the specular direction ${\hat{\textbf t}}$ at $s$ and then scattered into an off-specular direction ${\hat{\textbf u}}$ due to a perturbation of the normal at $s$ described probabilistically using $p$. This interpretation was explored further in our previous work [22].

We now introduce so-called global angles. By this, we mean that they are measured with respect to a static coordinate system rather than the local unit normal ${\hat{\textbf n}}$—see Fig. 6. In particular, let $\beta$ be the direction of ${\hat{\textbf n}}$ at $s$, let $\psi$ be the specular direction, and let $\gamma$ be the scattered direction in the global coordinate system, i.e.,

 

Fig. 6. Relation between the local and global angles; ${\vartheta _i} \lt 0$ and $\beta$, $\gamma$, $\psi$, ${\vartheta _r}$, ${\vartheta _o} \gt 0$.

Download Full Size | PDF

By inspecting Fig. 6, we see that (after some trigonometry)

Transforming ${\vartheta _o}$ to $\gamma$ and ${\vartheta _r}$ to $\psi$ in Eq. (23) yields

Suppose that we are interested in all light scattered by the surface in some global direction $\gamma$. In this case, we must integrate over the whole reflector, keeping in mind that the normal ${\hat{\textbf n}}$, and hence $\beta$, changes as we traverse the reflector. That is,

Notice that the left-hand side only depends on $\gamma$ so that we can define it as $h(\gamma)$ [${\rm W} \cdot {{\rm rad}^{- 1}}$]. It represents the scattered intensity distribution in direction $\gamma$. By changing the order of integration on the right-hand side, we get

This is a convolution integral, which is often denoted with an asterisk operator; Eq. (29), thus, has the shorthand notation $h(\gamma) = (p*g)(\gamma)$. We will occasionally use this notation, or simply $h = p*g$, referring to the convolution integral in Eq. (29). Note that this is a general relation in that it is independent of the light source, which is intuitive since the scattering event occurs at the reflector surface, irrespective of how the light arrives there.

We will discuss in detail how this model can be used to design specular reflectors while taking surface scattering into account in Section 3. However, the key is to prescribe the desired target distribution $h$ and the PDF $p$ describing the orientations of the microfacets. Next, a simple deconvolution procedure can be used to obtain the virtual specular intensity distribution $g$, which may be utilized to compute the reflector surface using well-established specular design procedures. We note that, due to surface scattering, the $g$ computed via deconvolution is not achieved in practice by the reflector manufactured using the scattering surface. Instead, the deconvolved $g$ is merely used as a tool to utilize well-established specular reflector design methods for designing reflectors with surfaces exhibiting scattering. Of course, a perfect mirror in the shape of the reflector computed using the deconvolved specular distribution would achieve this $g$.

2. Parallel Sources

Before discussing the practical considerations of using this theory for reflector design by showing some numerical examples, let us briefly derive the relations needed for specular reflector design, starting with parallel sources.

Recall Fig. 2 and Definition 1, where the parallel source was first introduced. In particular, consider a parallel source with ${\vartheta _s} = 0$ so that the rays travel along the positive $y$ axis. Furthermore, let $x$ be the coordinate running along the source, and let ${\ell _s}$ be the length of the source. The reflector is parametrized by ${\textbf r}(x) = (x,u(x{))^ \intercal}$, where $u(x) \gt 0$ for all $x \in [0,{\ell _s}]$ is a smooth height function. Thus, a unit tangent vector ${\hat{\boldsymbol \tau}}$ to the reflector is by definition given by ${\hat{\boldsymbol \tau}} = {{\textbf r}^\prime}(x)/| {{{\textbf r}^\prime}(x)} | = (1,{u^\prime}(x{))^ \intercal}/\sqrt {1 + {u^\prime}{{(x)}^2}}$, and the unit normal ${\hat{\textbf n}}$ pointing toward the light source is given by

The vectorial law of reflection may be expressed as ${\hat{\textbf t}} = {\hat{\textbf s}} - 2({\hat{\textbf s}} \cdot {\hat{\textbf n}}){\hat{\textbf n}}$ so that

Now, from the geometry near the intersection point $s$ between a source ray and the reflector surface (Fig. 7), we get

 

Fig. 7. Parallel illumination of a reflector traced by the curve $y = u(x)$; the relevant unit vectors and angles close to the intersection point $s$ are shown; ${\vartheta _i} \lt 0$ and $\psi$, ${\vartheta _r} \gt 0$.

Download Full Size | PDF

Then, by the definition of the dot product,

Returning to Eq. (32), we get

Finally, using Eq. (31), we find that

Solving for ${u^\prime}(x)$ gives the ODEs,

The sign of the root is related to the direction of the reflected ray. In particular, using the vectorial law of reflection, ${\hat{\textbf t}} = {\hat{\textbf s}} - 2({\hat{\textbf s}} \cdot {\hat{\textbf n}}){\hat{\textbf n}}$, we get ${\hat{\textbf t}} = (\pm \sin (\psi), - \cos (\psi {))^ \intercal}$ depending on the sign of the root. Since ${\hat{\textbf t}} = (\sin (\psi), - \cos (\psi {))^ \intercal}$, we use the positive root together with the tangent half-angle relation ([23], p. 127) to get the appropriate ODE,

To obtain the height function $u(x)$, we must solve this ODE with some initial value for the distance of the first point on the reflector from the source, $u(0) = {u_0}$, making it an initial value problem (IVP). Before discussing how to solve this IVP, note that $\psi$ is a map from a position on the source $x$ onto a direction $\psi$. Denoting this map $m(x)$, i.e., $\psi = m(x)$, we get the following IVP:

Note that a representation of the solution to this IVP is

Optical map. We will now show how to construct two types of solutions to the above IVP. In particular, suppose the optical map $m(x) = {m_ +}(x)$ is monotonically increasing on $x \in [0,{\ell _s}]$. Let $\forall\; x \in [0,{\ell _s}]$: $M(x) \gt 0$ [${\rm W} \cdot {{\rm m}^{- 1}}$] be the parallel source exitance, and recall that $\forall\; \psi \in [{\psi _1},{\psi _2}] \subseteq (- \pi /2,\pi /2)$: $g(\psi) \gt 0$ [${\rm W} \cdot {{\rm rad}^{- 1}}$] is the specular intensity. Since ${m_ +}(x)$ is monotonically increasing, conservation of flux gives

On the other hand, suppose the optical map $m(x) = {m_ -}(x)$ is monotonically decreasing on $x \in [0,{\ell _s}]$, then the conservation of flux gives

3. Point Sources

We now move on to point sources; recall Fig. 2 and Definition 2, where they were first introduced. In particular, consider a point source with ${x_s} = 0$, i.e., centered at the origin of an $x - y$ coordinate system. Let $\vartheta$ be the angle between the positive $y$ axis and a source ray—see Fig. 8. Let $[{\vartheta _1},{\vartheta _2}] \subseteq (- \pi /2,\pi /2)$ be the interval where the reflector, described using a smooth, continuous radius function $u(\vartheta) \gt 0$, exists. The reflector is, thus, parametrized by ${\textbf r}(\vartheta) = u(\vartheta){{\hat{\textbf e}}_r}$, where ${{\hat{\textbf e}}_r} = (- \sin (\vartheta),\cos (\vartheta {))^ \intercal}$ is the radial unit vector in our polar coordinate system. Thus, a unit tangent vector ${\hat{\boldsymbol \tau}}$ to ${\textbf r}(\vartheta)$ is given by ${\hat{\boldsymbol \tau}} = {{\textbf r}^\prime}(\vartheta)/| {{{\textbf r}^\prime}(\vartheta)} | = ({u^\prime}(\vartheta){{\hat{\textbf e}}_r} + u(\vartheta){{\hat{\textbf e}}_\vartheta})/\sqrt {u{{(\vartheta)}^2} + {u^\prime}{{(\vartheta)}^2}}$, where ${{\hat{\textbf e}}_\vartheta} = {\hat{\textbf e}}_r^\prime = (- \cos (\vartheta), - \sin (\vartheta {))^ \intercal}$ is the angular unit vector in our polar coordinate system. Whence, the unit normal ${\hat{\textbf n}}$ pointing toward the light source is given by

 

Fig. 8. Point source illumination of a reflector traced by the curve ${\textbf r}(\vartheta) = u(\vartheta){{\hat{\textbf e}}_r}$; the relevant unit vectors and angles close to the intersection point $s$ are shown; $\vartheta$, ${\vartheta _i} \lt 0$ and $\psi$, ${\vartheta _r} \gt 0$.

Download Full Size | PDF

Introducing $v(\vartheta): = \ln (u(\vartheta))$, we get

Note that all source rays ${\hat{\textbf s}}$ are given by ${\hat{\textbf s}} \equiv {{\hat{\textbf e}}_r}$ so that

From the geometry near the intersection of a source ray and the reflector surface (Fig. 8) and the definition of the dot product, we get, after some trigonometry (note that $\psi - \vartheta = {\vartheta _r} - {\vartheta _i}$ and recall that ${\vartheta _i} = - {\vartheta _r}$),

Since the law of reflection, ${\hat{\textbf t}} = {\hat{\textbf s}} - 2({\hat{\textbf s}} \cdot {\hat{\textbf n}}){\hat{\textbf n}}$, is valid for point sources too, so is Eq. (32), and we get

Finally, using Eq. (47), we find that

Solving for ${v^\prime}(x)$ gives the ODEs,

Note that, this time, $\psi$ is a map from the angle of a source ray, $\vartheta$, onto a specular direction $\psi$, i.e., $\psi = m(\vartheta)$. Comparing ${\hat{\textbf t}} = (\sin (\psi), - \cos (\psi {))^ \intercal}$ with ${\hat{\textbf t}} = {\hat{\textbf s}} - 2({\hat{\textbf s}} \cdot {\hat{\textbf n}}){\hat{\textbf n}} = - \cos (\vartheta - \psi){{\hat{\textbf e}}_r} \pm \sin (\vartheta - \psi){{\hat{\textbf e}}_\vartheta}$, we conclude that we must once again use the positive root. By using the tangent half-angle relation ([23], p. 127), we, thus, get

Using the initial condition $u({\vartheta _1}) = {u_0}$, we get the IVP,

After solving the IVP for $v(\vartheta)$, we obtain the radius function $u(\vartheta) = \exp (v(\vartheta))$, which fully determines the reflector.

Optical map. Similar to the parallel source case, we will now show how to construct two types of solutions to the above IVP. Suppose the optical map $m(\vartheta) = {m_ +}(\vartheta)$ is monotonically increasing on $\vartheta \in [{\vartheta _1},{\vartheta _2}] \subseteq (- \pi /2,\pi /2)$. Let $\forall\; \vartheta \in [{\vartheta _1},{\vartheta _2}]$: $f(\vartheta) \gt 0$ [${\rm W} \cdot {{\rm rad}^{- 1}}$] be the point source intensity, and recall that $\forall\; \psi \in [{\psi _1},{\psi _2}] \subseteq (- \pi /2,\pi /2)$: $g(\psi) \gt 0$ [${\rm W} \cdot {{\rm rad}^{- 1}}$] is the specular intensity. Since ${m_ +}(\vartheta)$ is monotonically increasing, conservation of flux gives

On the other hand, supposing the optical map $m(\vartheta) = {m_ -}(\vartheta)$ is monotonically decreasing on $\vartheta \in [{\vartheta _1},{\vartheta _2}]$, then conservation of flux gives

3. NUMERICAL EXAMPLES

This section shows how to use the equations we just derived. In particular, we will discuss the proposed workflow of an optical engineer and show how to solve the relevant differential equations for a selection of examples. For more examples and further discussion, please see our previous work [22].

A. Verification

To verify our model, the curves representing the reflectors were first computed numerically on a fixed grid (either in $x$ or $\vartheta$, depending on the source type) by solving the IVPs for $m$ and $u$. Next, raytracing was performed with rays sampled from the desired source distribution $M(x)$ or $f(\vartheta)$. In the custom raytracing procedure, a sampled ray strikes the reflector at some point $s$, where the normal ${\hat{\textbf n}}$ is determined using linear interpolation between known normals (i.e., at the fixed grid points where the reflector surface and normals were computed) using MATLAB’s interp1 routine. The law of reflection, ${\hat{\textbf t}} = {\hat{\textbf s}} - 2({\hat{\textbf s}} \cdot {\hat{\textbf n}}){\hat{\textbf n}}$, is subsequently applied to get the specular direction at $s$. Finally, the scattered ray ${\hat{\textbf u}}$ is computed by applying a rotation ${\textbf R}(\eta)$ with $\eta \in [{\eta _1},{\eta _2}] \subseteq (- \pi /2,\pi /2)$ sampled from the PDF $p(\eta)$, to ${\hat{\textbf n}}$ to get ${{\hat{\textbf n}}^m} = {\textbf R}(\eta){\hat{\textbf n}}$, and then the law of reflection, ${\hat{\textbf u}} = {\hat{\textbf s}} - 2({\hat{\textbf s}} \cdot {{\hat{\textbf n}}^m}){{\hat{\textbf n}}^m}$, is applied.

The spatial or angular domains are divided into $N$ equispaced regions, so-called “bins,” i.e., for some positive integer $j \lt N$: $[{x_{j - 1}},{x_j}]$ or $[{\vartheta _{j - 1}},{\vartheta _j}]$. The bins of the source, specular, and scattered rays are identified using MATLAB’s dsearchn nearest point search, and finally the ray counts are converted to an exitance $E$ or intensity $I$ using

Naturally, Eq. (59) is only used for the source rays from a parallel source since the specular and scattered rays all result in an intensity, which exists in angular space. For the specular and scattered rays with a parallel source, the following expression is used:

In other words, for a parallel source, $E$ approximates $M$ via Eq. (59), and $I$ approximates $g$ and $h$ via Eq. (61), while for a point source, $I$ approximates $f$, $g$, and $h$ via Eq. (60). More details are available in ([24], p. 34).

To quantitatively evaluate the accuracy of our predicted scattered distributions, we made use of the RMS error given by

B. Parallel Source

We will start with an example of reflector design using a parallel source. The first example is shown below, where

Example #1: Overlapping Gaussians:

Since the specular target distribution $g(\psi)$ is a sum of two Gaussians, it does not have a finite support. In order to clamp the support $[{\psi _1},{\psi _2}] \subseteq (- \pi /2,\pi /2)$, we fixed $\epsilon ={ 10^{- 3}}$ and found the points where $g(\psi) = \epsilon$. We then set the values of $g(\psi)$ to zero outside $[{\psi _1},{\psi _2}]$ and renormalized it such that energy was conserved. The Gaussian distribution $p(\eta)$ was given finite support $[{\eta _1},{\eta _2}] = [- \pi /2,\pi /2]$ in a similar manner, and it was then renormalized to unity on this domain, as required for it to be a PDF. The predicted scattered distribution $h(\gamma)$ was computed using MATLAB’s $\texttt{conv}$ function.

The distributions in this example are shown in Fig. 9. Note that all distributions are plotted on a range exceeding the previously mentioned domains. This is called “zero-padding,” which is necessary when performing convolution since the scattered distribution broadens so that it extends outside the support of $g$. When numerically computing the reflectors, we have chosen 256 sample points for our distributions. The reason for using numerical methods here is that, while we do have exact expressions for all the distributions in this problem, we want the method to be as general as possible.

 

Fig. 9. Distributions in Example #1; 256 sample points.

Download Full Size | PDF

To proceed, we now show the finite support $[{\psi _1},{\psi _2}]$ of $g(\psi)$ (found using $g(\psi) = \epsilon$), the specular mappings ${m_ \pm}$, and the computed reflectors in Fig. 10 for 256 sample points of the distributions. The specular mappings ${m_ \pm}$ were computed on a variable grid using MATLAB’s $\texttt{ode15s}$ solver with tolerances RelTol $={ 10^{- 12}}$ and AbsTol $={ 10^{- 14}}$, while the reflectors were computed using the integral representation of the solution, i.e., Eq. (40), using MATLAB’s integral command, on a fixed grid with 256 sample points.

 

Fig. 10. Finite domain of $g$, specular mappings ${m_ \pm}$, and associated reflectors ${{\textbf r}_ \pm}$ in Example #1; 256 sample points of $g$ and ${{\textbf r}_ \pm}$, initial condition ${u_ \pm}(0) = 1$.

Download Full Size | PDF

 

Fig. 11. Raytraced distributions in Example #1; 256 sample points, ${10^6}$ rays traced.

Download Full Size | PDF

We then raytraced these reflectors to validate our model of scattering. This procedure was outlined in Section 3.1, and the result is shown in Fig. 11 for the ${{\textbf r}_ -}$ reflector, computed using the monotonically decreasing optical map ${m_ -}$. Note the near-perfect ${-}1/2$ exponential convergence when the RMS error $\varepsilon (h,{h^*})$—recall Eq. (62)—is plotted against the number of rays traced, ${N_{\rm r}}$, as expected of Monte Carlo raytracing ([24], p. 9). The only problematic areas in the raytraced distributions are those close to ${\psi _1}$ and ${\psi _2}$ for ${g^*}$. This is smoothed out by the effect of scattering so that the scattered distribution is nearly ideal. As this manuscript is focused on scattering in freeform reflector design, we are content with this minor discrepancy.

C. Point Source

1. Example #2: Unknown Specular Distribution

Having verified our model, we will now consider an example problem using a point source—see the example below. This differs from the previous example in two crucial ways, in addition to the apparent difference in symmetry going from a parallel source to a point source. Namely, the exact $g$ is no longer known, as we only prescribed a desired scattered target distribution $h$. In addition, the scattering function $p$ is now given by a Lorentzian (also known as a Cauchy distribution),

Example #2: Unknown specular distribution:

The distributions are shown in Fig. 12, where we note the sharp cutoff of the tails of $p$ at $\eta = - \pi /2$ and $\eta = \pi /2$. The distributions underwent the same treatment as in the first example, i.e., clamping to finite domains and renormalization as appropriate.

 

Fig. 12. Distributions in Example #2; 256 sample points.

Download Full Size | PDF

 

Fig. 13. Finite domain of ${g_{{\rm dc}}}$, specular mappings ${m_ \pm}$, and associated reflectors ${{\textbf r}_ \pm}$ in Example #2; 256 sample points of ${g_{{\rm dc}}}$ and ${{\textbf r}_ \pm}$, initial condition ${u_ \pm}({\vartheta _1}) = 1$.

Download Full Size | PDF

 

Fig. 14. Raytraced distributions, $h$ and ${h_{{\rm rc}}}$ in Example #2; 256 sample points, ${10^6}$ rays traced.

Download Full Size | PDF

To compute the reflectors, we first determined a deconvolved specular target distribution ${g_{{\rm dc}}}$ using 10 iterations of MATLAB’s deconvlucy deconvolver, implementing the Richardson–Lucy deconvolution scheme. We then solved the IVPs for the optical mappings ${m_ \pm}$ as well as those for the reflector radius functions ${u_ \pm}$. This is shown in Fig. 13, together with the limits ${\psi _1}$ and ${\psi _2}$ (found via ${g_{{\rm dc}}}(\psi) = \epsilon$) used to solve the IVPs for ${m_ \pm}$—recall Eqs. (56) and (58). The IVPs for ${m_ \pm}$ and ${u_ \pm}$ were solved in the same way as in the previous example, i.e., using the ODE solver ode15s with improved tolerances and by solving the equivalent integral equation, Eq. (54), using MATLAB’s integral command, respectively.

We then raytraced the ${{\textbf r}_ +}$ reflector to obtain the data in Fig. 14, where we again observe the expected ${-}1/2$ slope of the RMS error. Looking at the specular distribution, we notice that the behavior close to the boundaries ${\psi _1}$ and ${\psi _2}$ is slightly worse than in the previous example. Given the fact that deconvolution is an ill-posed problem ([25], Ch. 1, Sec. V, p. 32), it is not obvious that the model would work well for this example. Indeed, by prescribing only the scattered distribution $h$, rather than the specular $g$, we are guaranteed neither existence nor uniqueness of ${g_{{\rm dc}}}$. Numerically, we can still apply MATLAB’s deconvlucy iterative deconvolver to find $a$ deconvolved distribution ${g_{{\rm dc}}}$, of course, but the existence of a true $g$ is not guaranteed. This is also why we compare the raytraced scattered distribution to ${h_{{\rm rc}}}: = p*{g_{{\rm dc}}}$, i.e., the “reconvolved” distribution, rather than the prescribed $h$. The difference between the two is shown at the bottom of Fig. 14, where it is possible to observe that ${h_{{\rm rc}}}$ extends further out when compared to $h$. This suggests certain restrictions on how close to an arbitrary target one can get within this framework, depending on the prescribed target distribution $h$ and surface scattering function $p$.

 

Fig. 15. Example #3; the effect of increased scattering on the radius of the reflectors; 256 sample points.

Download Full Size | PDF

2. Example #3: Effect of Scattering on Reflector

In order to get an appreciation of the effect scattering has on the shape of the reflector, we will now consider the problem in the box below. The setup is as follows: an isotropic point source $f$ is to be transformed into a broad Gaussian distribution $h$. We will first compute the specular reflectors that achieve the desired target $h$ by choosing $g \equiv h$ and subsequently introduce more scattering in the system. The surface scattering function $p$ was chosen to be a Gaussian with varying standard deviation ($\sigma$) depending on how rough the surface was. Note that the convolution between two Gaussians is a broader Gaussian ([25], p. 11). Thus, assuming $h$ is broad enough, we are guaranteed the existence of ${g_{{\rm dc}}}$ for each $\sigma$ (see Example #2 showing what can happen if we are not guaranteed the existence of ${g_{{\rm dc}}}$).

Example #3: Effect of scattering on reflector:

To compute the base reflectors ($\sigma = 0$) with initial condition ${u_ \pm}({\vartheta _1}) = 1$, the light distribution $h$ was used as the specular target distribution, i.e., $g \equiv h$, and the values of ${\psi _1}$, ${\psi _2}$ were determined as before. For each successive set of reflectors, the scattered target distribution was kept as $h$, but the changing $\sigma$ in $p$ resulted in different deconvolved specular target distributions ${g_{{\rm dc}}}$ for each value of $\sigma$. As before, these specular target distributions were used to compute the reflectors.

Next, the radius of each sample point on the reflectors was computed and compared to the radius of the base reflectors, thus producing the graphs in Fig. 15. We found that the differences in terms of radius are in the order of a few percent, depending on $\sigma$ and varying along the reflector. Note that the reflectors all start at the fixed point ${u_ \pm}({\vartheta _1}) = 1$ (initial condition) and, in this case, they return to a fixed point at ${u_ \pm}({\vartheta _2})$ as a consequence of the symmetry of the distributions. Since the problem of computing the reflectors is scale invariant, the values of a few percent difference in radius apply regardless of the unit of length in the system. The fact that all the computed reflectors achieved the desired target distribution (when taking scattering into account) was verified using raytracing (not shown).

This example demonstrates the added freedom in terms of designs that an approach such as ours can provide to an optical engineer since they can specify not only which target they would like to achieve but how they would like to achieve it. That is, they can consider the trade-off between added scattering and reduced glare versus, e.g., transfer efficiency, depending on the relevant design parameters of the problem.

4. CONCLUSION

We have presented a general method of computing two-dimensional scattering reflector surfaces. Specifically, by considering all the scattered light from a curve endowed with microfacets, a convolution integral between a PDF related to the stochastic orientation of the microfacets—and thus to surface roughness—and a virtual specular target distribution representing the reflected light from a smooth surface was reached. Deconvolution, thus, yields an approximation of the specular target distribution after prescribing a desired scattered target light distribution and the PDF. In this way, the problem is reduced to a well-understood specular inverse problem. In its current form, which we view as a proof of principle as we work our way toward more realistic unconstrained three-dimensional scattering, the main strength of our approach is, thus, that it can be performed as a pre-processing step before using already mature specular design workflows.

The model was verified by comparing the predicted light distribution from the convolution integral with one obtained using a custom raytracer, which directly implemented the microfacets used as a basis for our model. This verification step showed that our model was consistent for parallel and point sources. We also found that scattering resulted in a change in reflector radius of a few percent versus the specular case, depending on the amount of surface scattering present in the system. Due to the ill-posedness of deconvolution, the proposed algorithm is most effective for smooth scattered target distributions.

This work makes it easy to think of many exciting extensions. Some examples include treating refractive media (i.e., lenses), three-dimensional freeform reflectors, and allowing the PDF related to the roughness of the surface to change as a function of position or to take into account effects such as shadowing and masking, i.e., the idea that adjacent microfacet may obstruct the flux incident or scattered from a given microfacet. We are actively working toward a more general three-dimensional model of surface scattering with the ultimate goal of designing three-dimensional freeform reflectors.