{'fraction_non_alphanumeric': 0.04382623076695904, 'fraction_numerical': 0.009361802663831546, 'mean_word_length': 4.441426964486582, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Two-Dimensional Freeform Reflector Design with a Scattering Surface

We combine two-dimensional freeform reflector design with a scattering surface modelled using microfacets, i.e., small specular surfaces representing surface roughness. The model results in a convolution integral for the scattered light intensity distribution, which yields an inverse specular problem after deconvolution. Thus, the shape of a reflector with a scattering surface may be computed using deconvolution, followed by solving the typical inverse problem of specular reflector design.


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][2][3][4][5]. A common situation involves shaping light from a point source since this is the approximate distribution of a light-emitting diode (LED). Romijn extensively studied shaping point source illumination using reflectors and lenses by constructing numerical solutions to the so-called generalised 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 a strong demand for aesthetic lighting in terms of colour and light shaping. Diffusive media or scattering surfaces can be utilised to homogenise 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, 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. 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 utilised or detrimental effects due to surface defects to be pre-emptively mitigated. Literature concerning the unification of scattering and freeform optics is minimal. The best reference we have found is Lin et al., who designed a lens with a freeform scattering inner surface and a spherical outer one [10]. The freeform surface was designed by an iterative optimisation procedure whereby a freeform surface represented by Bézier curves was first computed and then modified iteratively to take into account the difference between the prescribed target distribution and a raytraced one.
In contrast, we shall focus on computing two-dimensional freeform reflectors with a scattering surface using a direct inverse method. Specifically, we shall 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 dictating the microfacets 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 roughness of the surface can be estimated a priori, deconvolution can be utilised to compute a virtual specular light distribution, which can then 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 -Sec. 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.

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, and 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 a cone of 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.

Geometry
The two-dimensional geometry of scattering used throughout this manuscript is outlined in Fig. 1. Here, we see an incident ray along directionŝ striking the surface at a point P, resulting in a specular ray alongt for perfect specular reflection or a scattered ray alongû. The directionsŝ,t andû are given by angles ϑ i ∈ (−π/2, π/2), ϑ r = −ϑ i and ϑ o ∈ (−π/2, π/2), respectively, i.e.,

Scattering Model
Having introduced the geometry, let us now consider how the model is constructed. We shall first introduce the concept of zero-étendue sources, followed by the bidirectional reflectance distribution function (BRDF) and, finally, the microfacets used to model surface roughness.

Zero-étendue Sources
The concept of étendue is easiest to conceptualise 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(ϑ) of the ray at x in direction ϑ, where n is the refractive index. We shall always take n as unity, so that p ≡ sin(ϑ).
The étendue in our two-dimensional phase space is given by, cf. [11, Eq. (5.50), p. 89], where X and Θ are the spatial and angular domains of integration, respectively. We shall consider two-dimensional parallel-ray sources (referred to as parallel sources henceforth) and point sources in this manuscript. In the case of parallel sources, the rays are emitted in a fixed direction ϑ s , modulated along the spatial extent s of the source, whilst in the case of point sources, the source is located at a fixed position x s , where it emits modulated light in all directions ϑ. The sources are illustrated in Fig. 2, where the point source has been restricted to emitting rays in the interval ϑ ∈ [ϑ 1 , ϑ 2 ] ⊆ (−π/2, π/2). The radiances L [W · rad −1 · m −1 ] of these sources are summarised in the definitions below.
Definition 1. The radiance of a two-dimensional parallel source emitting rays in direction ϑ s is given by where M (x) [W · m −1 ] is the exitance of the source at position x and δ(ϑ − ϑ s ) [rad −1 ] is the Dirac delta function.
Definition 2. The radiance of a two-dimensional point source at position x s is given by where I(ϑ) [W · rad −1 ] is the intensity of the source in direction ϑ and δ(x − x s ) [m −1 ] is the Dirac delta function.
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 shall consider in this manuscript.

Bidirectional Reflectance Distribution Function
The bidirectional reflectance distribution function (BRDF) B(ϑ i , ϑ o ) [rad −1 ] was first formulated by Nicodemus in 1965 [12], and it describes the way light from direction ϑ i is scattered from the surface into direction ϑ o [13] (recall Fig. 1). To see how, let L i (ϑ i ) [W · rad −1 · m −1 ] be the incoming radiance from direction ϑ i , and let It can be shown that any physical BRDF must satisfy the following properties [14]: • Positivity: which must hold for obvious reasons -both L i and L o are nonnegative in Eq. (5), so naturally, this applies to the BRDF as well.
• Helmholtz reciprocity: which encapsulates the intuitive fact that the flux in direction ϑ o from a source at ϑ i is the same as the flux in direction ϑ i from the same source at ϑ o , assuming the rest of the system remains the same.
• Energy conservation: ∀ϑ i ∈ (−π/2, π/2) : which represents conservation of flux, implying that all the incoming light must be scattered in some direction. Note that this could be less than unity if absorbance or other losses were considered, but these effects are outside this work's scope.
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 ϑ o ≡ ϑ r = −ϑ i since the angles originate from the unit surface normaln at the point of intersection. This fact leads to the following definition for the BRDF of specular reflection.
Definition 3. The bidirectional reflectance distribution function (BRDF) of specular reflection is given by where δ is the Dirac delta function.
This BRDF fulfils: • Positivity by the properties of the Dirac delta function.
• Helmholtz reciprocity, Eq. (7): since the Dirac delta function is nonzero only when ϑ o = −ϑ i , and cosine is even.
Additionally, the scattering equation, Eq. (5), yields Note that we are keeping the cosine terms in this relation for later, but note that we could write That is, the radiance in the outgoing direction ϑ o is exactly the radiance from the incoming direction ϑ i , which is what we expect for specular reflection.

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 probability density function, a realistic model for surface roughness may be constructed. This concept has been used extensively in computer-generated imagery (CGI) [15,Ch. 8.4]. Of particular interest is the microfacet model developed by Torrance and Sparrow in 1967 [16], and later applied to CGI by Cook and Torrance in 1982 [17]. We have opted to develop a simpler model inspired by the one introduced by Torrance and Sparrow (in particular, we will 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 whilst remaining realistic.
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,n, aligns with the y-axis; this yields the view in Fig. 4 and defines our coordinate system with its origin at P. Finally, consider tilting this locally flat reflector piece by some η ∈ (−π/2, π/2), whilst keeping the angle of the incident light constant with respect to the original unit normaln -see Fig. 5, where ϑ m i and ϑ m o are the incident and outgoing directions with respect to the microfacet normaln m , respectively. Mathematically, we define a rotation matrix to construct the microfacet normaln m := R(η)n. Since the microfacet is specular, the law of reflection yields ϑ m o ≡ ϑ m r = −ϑ m i . 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 two times the perturbation of the normal. Inserting this expression for ϑ 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 η with respect to the macroscopic surface normal.
Distributed orientations To use microfacets to describe scattering from a surface, rather than fixing η, we shall 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 probability density functionp(η) with support supp(p) = [η 1 , η 2 ] ⊆ (−π/2, π/2), i.e., η2 η1p (η) dη = 1. Note that whilst the support ofp can formally extend to (−π/2, π/2), we will typically consider situations where the significant contributions are concentrated in a much more narrow region centred around η = 0. This shape is consistent with measured BRDFs from machined surfaces [18,Ch. 4].
To get the outgoing radiance, we now multiply the right-hand side of Eq. (15) byp(η) and integrate over all valid microfacet orientations, leading to the below definition. Definition 4. Let η ∈ [η 1 , η 2 ] ⊆ (−π/2, π/2) be the microfacet orientation with respect to the macroscopic normal n, sampled from the probability density functionp(η). Let L o be the outgoing (i.e., scattered) radiance in direction ϑ o due to an incoming radiance L i from direction ϑ i . Then: By using η ≡ (ϑ i + ϑ o )/2 -recall Eq. (14) -and changing the integration variable to ϑ i , we get Since ϑ o ∈ (−π/2, π/2), and supp(p) = [η 1 , η 2 ] ⊆ (−π/2, π/2), we can rewrite this integral as Note that the above expression does not consider whether a ray can physically reach ϑ o from ϑ i . As an extreme example, consider near-grazing incident light, where any nonzero perturbation η 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 α : where which is the so-called microfacet scattering equation for microfacets with orientations described by the probability density function p, via p(2η) =p(η)/2.

Freeform Reflector Design
Thus far, we have derived the scattered radiance [W · rad −1 · 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 shall show how this situation yields a convolution integral for the outgoing intensity [W · 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 [8, Chs. 1 − 6].

Intensities
Let s ∈ [0, ], parametrising a curve of length , be the arc length. The point-wise microfacet scattering equation, Eq. (21), then becomes where ϑ o and ϑ i are measured from the unit normaln of the reflector at s, and L i (s, ϑ i ), L o (s, ϑ o ) are the incident and outgoing radiances at s, respectively. Note that we have assumed that the scattering is independent of the location along the reflector since p does not depend on s, i.e., the surface is assumed to be isotropic. For a perfectly smooth reflector, an incident ray is specularly reflected from direction ϑ i into direction ϑ r = −ϑ i . Defining a virtual specularly reflected radiance L r (s, ϑ r ) := L i (s, −ϑ 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 t at s and then scattered into an off-specular directionû due to a perturbation of the normal at s described probabilistically using p. This interpretation was explored further in our previous work [19]. 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 normaln -see Fig. 6. In particular, let β be the direction ofn at s, let ψ be the specular direction, and let γ be the scattered direction in the global coordinate system, i.e., n = sin(β), − cos(β) ,t = sin(ψ), − cos(ψ) ,û = sin(γ), − cos(γ) . (24) By inspecting Fig. 6, we see that (after some trigonometry) Transforming ϑ o to γ and ϑ r to ψ in Eq. (23) yields where we have used the fact that p has finite support [η 1 , η 2 ] ⊆ (−π/2, π/2) to extend the integration limits to the full circle. Suppose that we are interested in all light scattered by the surface in some global direction γ. In this case, we must integrate over the whole reflector, keeping in mind that the normaln, and hence β, changes as we traverse the reflector. That is, Notice that the left-hand side only depends on γ, so that we can define it as h(γ) [W · rad −1 ]. It represents the scattered intensity distribution in direction γ. By changing the order of integration on the right-hand side, we get dψ, where g(ψ) also has units [W · rad −1 ] and represents a virtual specular intensity distribution in direction ψ. In total, we thus have that the outgoing intensity h in some direction γ is given by This is a convolution integral, which is often denoted as an asterisk operator; Eq. (29) thus has the shorthand notation h(γ) = (p * g)(γ). 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, independently of how the light arrives there. We shall discuss in detail how this model can be used to design specular reflectors whilst taking surface scattering into account in Sec. 3. However, the key is to prescribe the desired target distribution h and the PDF p describing the orientations of the microfacets. Then, a simple deconvolution procedure can be used to obtain the virtual specular intensity distribution g, which is then utilised to compute the reflector surface using well-established specular design procedures. We note that due to surface scattering, the g we computed by 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 utilise 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.

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 ϑ s = 0 so that the rays travel along the positive y-axis. Furthermore, let x be the coordinate running along the source, and let s be the length of the source. The reflector is parametrised by r(x) = x, u(x) , where u(x) > 0 for all x ∈ [0, s ] is a smooth height function. Thus, a unit tangent vectorτ to the reflector is by definition given byτ = r (x)/|r (x)| = 1, u (x) / 1 + u (x) 2 , and the unit normaln pointing towards the light source is given bŷ where the rotation matrix R was introduced in Eq. (13). Recall that all source raysŝ are given byŝ ≡ê y = (0, 1) , so that The vectorial law of reflection may be expressed ast =ŝ − 2(ŝ ·n)n, so that Now, from the geometry near the intersection point s between a source ray and the reflector surface (Fig. 7), we get Then, by the definition of the dot product, −ŝ ·t = cos(2ϑ r ) = cos(ψ).
Finally, using Eq. (31), we find that − cos(ψ) = 1 − 2 Solving for u (x) gives the ordinary differential equations (ODEs) The sign of the root is related to the direction of the reflected ray. In particular, using the vectorial law of reflection,t =ŝ − 2(ŝ ·n)n, we gett = (± sin(ψ), − cos(ψ)) depending on the sign of the root. Sincê t = (sin(ψ), − cos(ψ)) , we use the positive root together with the tangent half-angle relation [20, 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 ψ is a map from a position on the source x onto a direction ψ. Denoting this map m(x), i.e., ψ = 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( be the parallel source exitance, and recall that ∀ψ ∈ [ψ 1 , ψ 2 ] ⊆ (−π/2, π/2): g(ψ) > 0 [W · rad −1 ] is the specular intensity. Since m + (x) is monotonically increasing, conservation of flux gives ∀x ∈ [0, s ] : which can be differentiated with respect to x to give the equivalent initial value problem: On the other hand, suppose the optical map m(x) = m − (x) is monotonically decreasing on x ∈ [0, s ], then conservation of flux gives ∀x ∈ [0, s ] : which can be differentiated with respect to x to give the equivalent IVP: (44)

Numerical Examples
This section shows how to use the equations we just derived. In particular, we shall 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 [19].

Verification
We have written a custom raytracer to verify our model, which implements the microfacets introduced when deriving the scattering model. First, the reflector surface was numerically computed on a fixed grid (either on x or ϑ, depending on the source type) by solving the IVPs for m and u. Then, the raytracing starts with rays sampled from the desired source distribution M (x) or f (ϑ). A sampled ray strikes the reflector at some point s, where the normaln 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, t =ŝ − 2(ŝ ·n)n, is subsequently applied to get the specular ray at s. Next, the scattered rayû is computed by applying a rotation R(η) with η ∈ [η 1 , η 2 ] ⊆ (−π/2, π/2) sampled from the PDF p(η), ton to getn m = R(η)n, and then the law of reflection,û =ŝ − 2(ŝ ·n m )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 < N : 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 or ∀j ∈ [1, N ] : for illumination by parallel and point sources, respectively. Here, the integral represents the total flux of the source, and Pr(x j−1 ≤ x < x j ) or Pr(ϑ j−1 ≤ ϑ < ϑ j ) is the number of rays in the jth bin divided by the total number of rays traced, i.e., the probability of falling in the jth bin, and ∆x or ∆ϑ is the size of the collection bins. 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), whilst for a point source, I approximates f , g and h via Eq. (60). More details are available in [21, p. 34].
To quantitatively evaluate the accuracy of our predicted scattered distributions, we shall use the root mean square (RMS) error given by for N collection bins of the raytraced scattered distribution, h * , and the scattered distribution h. We always denote raytraced distributions with an upper-index asterisk ( * ). In addition, the deconvolved specular distribution is denoted g dc , whilst h rc := g dc * p (recall Eq. (29)) represents the "reconvolved" distribution.

Parallel Source
We shall start with an example of reflector design using a parallel source. The first example is shown in the box below, where represents the Gaussian distribution centred at µ with standard deviation σ. The distributions M , g, and p have finite support given by [0, s ], [ψ 1 , ψ 2 ], and [η 1 , η 2 ], respectively. That is, we set the values to exactly zero outside these domains and renormalise as appropriate. Since we first wish to verify our predicted scattered distribution h, we shall prescribe the known specular distribution g and use this as our target function when designing the reflectors. Thus, the verification step here is to check that the raytraced scattered distribution h * approaches the predicted scattered distribution h from the convolution integral in Eq. (29).  Since the specular target distribution g(ψ) is a sum of two Gaussians, it does not have a finite support. In order to clamp the support [ψ 1 , ψ 2 ] ⊆ (−π/2, π/2), we fixed = 10 −3 and found the points where g(ψ) = . We then set the values of g(ψ) to zero outside [ψ 1 , ψ 2 ] and renormalised it such that energy was conserved. The Gaussian distribution p(η) was given finite support [η 1 , η 2 ] = [−π/2, π/2] in a similar manner, and it was then renormalised to unity on this domain, as required for it to be a probability density function. The predicted scattered distribution h(γ) was computed using Matlab's 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, whilst we do have exact expressions for all the distributions for this problem, we want the method to be as general as possible. To proceed, we now show the finite support [ψ 1 , ψ 2 ] of g(ψ) (found using g(ψ) = ), the specular mappings m ± , and the computed reflectors in Fig. 10, for 256 sample points of the distributions. The specular mappings m ± were computed on a variable grid using Matlab's ode15s solver with tolerances RelTol = 10 −12 and AbsTol = 10 −14 , whilst the reflectors were computed using the integral representation of the IVP, i.e., Eq. (40), using Matlab's integral command, on a fixed grid with 256 sample points.
We then raytraced these reflectors to validate our model of scattering. This procedure was outlined in Sec. 3.1 and the result is shown in Fig. 11 for the r − reflector, computed using the monotonically decreasing optical map m − . Note the near-perfect −1/2 exponential convergence when the RMS error ε(h, h * ) -recall Eq. (62) -is plotted against the number of rays traced, N r , as expected from Monte Carlo raytracing [21, p. 9]. The only problematic areas in the raytraced distributions are those close to ψ 1 and ψ 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 perfectly content with this minor discrepancy.

Point Source
Having verified our model, we shall now consider an example problem using a point source -see the box 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 prescribe a desired scattered target distribution h. In addition, the scattering function p is now given by a Lorentzian (also known as a Cauchy distribution) p(η) = L(η; σ) := 1 πσ with a full width at half maximum (FWHM) of 2σ. This change in p is significant for two reasons: machined mirrors often exhibit this type of BRDF [18,Ch. 4], and the tails fall to zero at a significantly lower rate compared to the Gaussian we used in our previous example.  The distributions are shown in Fig. 12, were we note the sharp cutoff of the tails of p at η = −π/2 and η = π/2. The distributions underwent the same treatment as in the first example, i.e., clamping to finite domains and renormalisation as appropriate. To compute the reflectors, we first determined a deconvolved specular target distribution g dc using ten iterations of Matlab's dconvlucy deconvolver, implementing the Richardson-Lucy deconvolution scheme. We then solved the IVPs for the optical mappings m ± as well as those for the reflector radius functions u ± . This is shown in Fig. 13, together with the limits ψ 1 and ψ 2 (found via g(ψ) = ) used to solve the IVPs for m ± -recall Eqs. (56) and (58). The IVPs for m ± and u ± 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. Figure 13: The finite domain of g dc , specular mappings m ± and associated reflectors r ± in Example #2; 256 sample points of g dc and r ± .
We then raytraced the r + reflector to obtain the data in Fig. 14, where we again observe the expected −1/2 trend with the number of rays traced. Looking at the specular distribution, we notice that the behaviour close to the boundaries ψ 1 and ψ 2 is slightly worse than in the previous example. Given the fact that deconvolution is an ill-posed problem [22, 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 not guaranteed existence or uniqueness of g dc . Numerically, we can still apply Matlab's deconvolucy iterative deconvolver to find a deconvolved distribution g dc , of course, but the existence of a true g is not a given. This is also why we compare the raytraced scattered distribution to h rc := p * g dc , i.e., the 'reconvolved' distribution, rather than the prescribed h. The difference between the two is shown at the bottom of

Conclusions
We have presented a general method of computing two-dimensional (and by extension rotationally and cylindrically symmetric) scattering reflector surfaces. Specifically, by considering all the scattered light from a curve endowed with microfacets, a convolution integral between a probability density function (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 virtual 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, 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, with reflectors computed using the deconvolved virtual specular target distribution. This verification step showed that our model was consistent for parallel and point sources. Due to the ill-posedness of deconvolution, the proposed algorithm is most effective for smooth scattered target distributions, and discontinuous scattered target distributions cannot be deconvolved since no sharper feature can be created.
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. In the case of refractive media, one could, for instance, use our model to correct for surface scattering due to machining imperfections. Allowing the PDF to change along the two-dimensional reflector results in a type of Fredholm integral equation. We are actively investigating every aforementioned extension.
Funding: This work was partially supported by the Dutch Research Council (Dutch: Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)) through grant P15-36. Disclosures: The authors declare no conflicts of interest. Data availability: Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon request.