Issue 
J. Eur. Opt. SocietyRapid Publ.
Volume 19, Number 1, 2023
EOSAM 2021



Article Number  18  
Number of page(s)  22  
DOI  https://doi.org/10.1051/jeos/2023014  
Published online  28 April 2023 
Research Article
Modelling surface light scattering for inverse twodimensional reflector design
^{1}
Department of Mathematics and Computer Science, Eindhoven University of Technology, PO Box 513, 5600 MB Eindhoven, The Netherlands
^{2}
Signify Research, High Tech Campus 7, 5656 AE Eindhoven, The Netherlands
^{*} Corresponding author: v.c.e.kronberg@tue.nl
Received:
27
January
2022
Accepted:
15
March
2023
We present a novel approach of modelling surface light scattering in the context of twodimensional reflector design, relying on energy conservation and optimal transport theory. For isotropic scattering in cylindrically or rotationally symmetric systems with inplane scattering, the scattered light distribution can be expressed as a convolution between a scattering function, which characterises the optical properties of the surface, and a specular light distribution. Deconvolving this expression allows for traditional specular reflector design procedures to be used, whilst accounting for scattering. This approach thus constitutes solving the inverse problem of light scattering, allowing for direct computation of the reflector surface, without the need for design iterations.
Key words: Surface scattering / Reflector design / Illumination optics / Inverse methods
© The Author(s), published by EDP Sciences, 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The engineering field of optical design is concerned with constructing an optical system that produces a desired output light distribution from a given source light distribution. This is often done within the confines of the socalled geometrical optics (GO) approximation, where light propagation is modelled using light rays – lines collinear with the Poynting vector. Within the GO approximation, phenomena such as diffraction and interference are typically not accounted for ([1], p. 159]). Contemporary optical designs sometimes incorporate scattering elements, but due to the absence of such phenomena in the GO approximation, their inclusion is typically realised using trialanderror in the design phase. This places restrictions on what kind of optical systems one can design, as scattering is difficult to account for in the design process. In this paper, we make a first effort towards constructing a surface light scattering model suitable for inverse reflector design.
There are a myriad of approaches one could consider to include scattering. The most fundamental, and in some ways most natural, approach is that of solving Maxwell’s equations directly. This would constitute a substantial departure in terms of strategy, and would in theory work for most realistic systems, but it is often impractical due to spiralling complexity. As such, a number of approximations based on Maxwell’s equations have been formulated which are applicable to problems within several regimes of parameter values, such as surface roughness or incident angle. Some highlights include the rigorous vector perturbation theory published by Lord Rayleigh in 1907 [2], and later expanded by Rice (1951) [3], and the rather different Kirchhoff approach, based on random phase variations due to microtopographic surface features, most commonly attributed to Beckmann and Spizzichino (1963) [4]. These approaches are valid in different regimes. In particular, Rayleigh–Rice vector perturbation theory agrees well with experimental measurements of wideangle scattering (up to approximately 50° of the polar angle of detector/source) for scattering from optically smooth surfaces. Here, “optically smooth” refers to the root mean square (RMS) surface roughness σ_{s} divided by the wavelength λ being much less than unity, i.e., σ_{s}/λ ≪ 1 ([5], p. 49). The Beckmann–Kirchhoff theory, on the other hand, is valid for rougher surfaces, but due to a moderateangle assumption as part of its derivation, it is not suitable for use with wide scattering angles and/or large angles of incidence.
There have been numerous developments since these early theories were formulated. Here, we highlight the work of Church, who published multiple papers during the 1970s on Rayleigh–Rice theory in the context of surface scattering from optically smooth surfaces [6]. According to Harvey, his contributions were instrumental in shaping the applied optics community at the time ([5], p. 50). The last work from this era we want to highlight is that of Harvey and Shack from 1976, where they developed a linear systems formulation of surface scattering based on a surface transfer function [7]. This approach allows for the use of the Fourier transform of the surface transfer function to compute a scattered radiance function “closely related to the bidirectional reflectance distribution function (BRDF)” ([5], p. 50). This was later extended to the generalised HarveyShack surface scattering theory, which is able to treat arbitrarily rough surfaces with arbitrary incident and scattering angles [8].
Rather than solving Maxwell’s equations, with or without approximations, one could alter the GO rays in some manner such that they can be used to compute scattering phenomena. One example of such an approach, that still retains many of the computational benefits of GO, is to modify the GO rays such that they carry information regarding the phase of the light, which can form the basis of diffraction calculations. A good overview of this approach can be found in McNamara [9]. Whilst computationally efficient, such an approach would still require substantial modifications to contemporary reflector design procedures.
In our approach, we remain in the domain of traditional GO, and as an alternative to carrying phase information, we propose a surface scattering model inspired by optimal transport theory [10], which leads to a convolution integral for cylindrically and rotationally symmetric problems with isotropic inplane scattering. This convolution integral yields the scattered light distribution, given a scattering function and a specular target distribution, where the former characterises the optical properties of the surface. This is a forward problem, but it can also be cast in terms of an inverse problem – given a desired target distribution and a scattering function, one may perform deconvolution to find an appropriate intermediate specular target distribution, which can in turn be used to design the optical element using traditional specular design methods. Our approach thus allows scattering to be taken into account in a simple preprocessing step, allowing for direct computation of the reflector surface without the need for design iterations or expensive forward raytracing calculations. This is in contrast with the established techniques of designing reflectors with a scattering surface, constituting forward scattering calculations and manual design iterations, i.e., designing a reflector and raytracing with a suitable Bidirectional Reflectance Distribution Function (BRDF) using software like LightTools, Zemax or Code V, followed by making manual adjustments to the reflector surface in an attempt to account for potential discrepancies between the resulting and desired scattered light distributions. The advantage of our approach is thus that we may separate the scattering calculations from the reflector design step, allowing us to greatly benefit from the maturity of specular reflector design procedures, whilst still accounting for scattering.
The computation of specular reflectors has a long history that is largely outside of the scope of this paper. A recent highlight is the development of methods based around solving the Monge–Ampère equation [11, 12]. A comprehensive literature review of inverse illumination optics may be found in [11] of Chapter 3. We shall restrict our attention to the twodimensional procedure developed by Maes ([13], Chap. 3) to compute cylindrically and rotationally symmetric reflectors. To the best of our knowledge, directly computing reflectors with scattering surfaces in illumination optics is very rare; 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 [14]. 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.
1.1 Summary of main contributions
In words, the main contribution of our approach is that it contains a decoupling between the inverse scattering problem and the inverse reflector computation. This is very powerful, since it allows one to take scattering into account in a preprocessing step, such as in the algorithm below.
Prerequisites
For our method, we require:

A source light distribution [lm/rad] (denoted f herein).

A prescribed target light distribution [lm/rad] (denoted h herein).

A surface scattering function [1/rad] (denoted p herein; related to the BRDF).
Algorithm

Deconvolve h and p to obtain a specular intensity distribution [lm/rad] (denoted g_{dc} herein).

Compute the specular reflector which transforms the source distribution f into g_{dc}.
The effect of each step is shown in Figure 1, with a point source at the origin. We note that the deconvolution requires some care, as a light distribution is necessarily nonnegative. For this reason, we used Richardson–Lucy deconvolution, which guarantees nonnegativity of the solution. Computing the specular reflector essentially reduces to solving two coupled initial value problems.
Figure 1 The inverse problem of reflector design with a scattering surface; the first step (deconvolution) yields the specular target distribution, which in turn can be used to compute the reflector by solving the inverse specular problem. The data is taken from Example #1 in Section 4. 
1.2 Structure of manuscript
The structure of this paper is as follows. Cylindrically symmetric problems with inplane scattering are covered first, together with a few words about deconvolution in Section 2.1, followed by rotationally symmetric problems with inplane scattering in Section 2.2. Next, twodimensional specular reflector design is briefly discussed in Section 3, followed by some results in the form of reflectors and raytraced distributions for validation of both the cylindrically and rotationally symmetric systems in Section 4. Finally, conclusions with some proposals for expansions of the model as well as suggestions for how to relate it to the BRDF are presented in Section 5.
2 Model
2.1 Cylindrically symmetric problems
Starting with cylindrically symmetric problems, consider the situation depicted in Figure 2. This specular problem can be viewed as a crosssection of a translationally invariant problem, such as an extruded optical element, with a linesource along the suppressed zaxis. For such a system, the specular problem may be analysed in twodimensions ([10], Chap. 3). Thus, the intensities and reflector surfaces are independent of z and we may study a crosssection in the plane of incidence, which is spanned by the source and reflected rays, along and , respectively, and which contains the unit normal . The hat () indicates unit vectors throughout this manuscript, and positive angles are by convention counterclockwise. The reflector is parametrised by , where u(φ) > 0 is at least twice continuously differentiable and is the radial unit vector in polar coordinates. The angle φ is measured counterclockwise from the positive axis, and it fully characterises the source ray along , emitted form the line source at the origin . The source ray along intersects the reflector at some point , where the unit normal of the reflector is given by . We take the convention , i.e., the normal is chosen directed towards the light source. From the specular law of reflection (LoR), we get an expression for the reflected direction , i.e.,(1)where we denote the angle between the positive axis and by ψ.
Figure 2 Reflectors exhibiting specular reflection (a), and light scattering (b); linesource along suppressed zaxis passing through . 
To introduce scattering, consider the situation depicted in Figure 2b. Inherent in this description is that we assume the scattering is limited to the plane of incidence, so that we can again study a crosssection of the translationally invariant problem. This assumption is made in order to preserve the twodimensional nature of the specular problem also when considering scattering. We postulate that it is valid for a situation where an extruded reflector is illuminated by a linesource parallel to the extrusion direction. This is because the process of extrusion leads to grooves along the extrusion direction, so that the crosssectional plane of symmetry is preserved. The validity of this claim has not been investigated further; it is thus an assumption made in the derivation of our twodimensional surface scattering model. Here, the source ray along gets mapped to a scattered ray along , where γ is measured counterclockwise from the positive xaxis. The scattered direction can be described as a rotation of by a stochastic parameter around the axis parallel to the zaxis passing through , i.e.,(2)
The stochastic parameter α is related to the scattering characteristics of the surface. We note that α depends on ψ, both in the sense that it will almost certainly have a different stochastic value for a given ψ – in fact, since α is sampled from a probability distribution, it has multiple values for all ψ — and in the sense that the probability distribution from which it is sampled may be different for different values of ψ. We shall return to the meaning of this, both mathematically and physically, in Section 2.1.2.
2.1.1 Mappings
To express reflection and scattering in terms of angles, let us introduce two mappings which give the reflected and scattered directions, i.e.,(3)where the former is a reformulation of the law of reflection, and the latter represents the scattering. This might seem superfluous, but it will simplify the discussion later, especially in the general threedimensional case, which we intend to treat in a future publication. In addition to these maps, we require their inverses to exist,(4)
Finally, we define a mapping yielding the scattering angle α, for fixed ψ and γ, i.e.,(5)
For a schematic summary, see Figure 3.
Figure 3 Relations between unit vectors and angles; all symbols are defined in the text. 
Explicitly, the maps are(6)where the first two relations follow from – see Figure 2a – and the LoR – Eq. (1). The existence of inverse mappings is not a priori guaranteed for all situations, but we shall restrict our attention to those where they do exist.
2.1.2 Energy balances
Having presented the mappings for the angles, we introduce the intensity distributions (illuminance) [lm/rad]:

source intensity distribution f(φ) > 0, φ ∈ [φ_{1}, φ_{2}],

intermediate specular intensity distribution g(ψ) > 0, ψ ∈ [ψ_{1}, ψ_{2}],

diffuse target intensity distribution h(γ) > 0, γ ∈ [γ_{1}, γ_{2}],
Consider next the relationship between ψ and γ for a fixed ψ = Ψ. Suppose we have a perfect specular reflector (i.e., a mirror). Then, α always vanishes, such that ψ ≡ γ and for fixed ψ = Ψ, we simply get a fixed γ = Γ. This is depicted schematically in Figure 4a. Suppose instead that we have nonzero scattering, and let the scattering vary depending on the incident angle. This yields a situation like the one depicted schematically in Figure 4b. Here, fixing ψ = Ψ and tracking where all the light emerges, we see that it falls within the interval .
Figure 4 Schematic maps for specular reflectors (a), and diffuse reflectors (b). 
Motivated by the above observation and the more fundamental concept of optimal transport theory, and in particular Monge–Kantorovich problems ([10], Chap. 1), let us introduce the density ρ(ψ, γ) > 0, ψ ∈ [ψ_{1}, ψ_{2}], γ ∈ [γ_{1}, γ_{2}], with properties(8a) (8b)that is, given ψ, integrating ρ(ψ, γ) over all scattered angles γ gives the intermediate specular intensity g(ψ) in the direction ψ. Stated otherwise, ρ(ψ, γ) determines how g(ψ) is spread over the range of scattered angles γ. Vice versa, given γ, integrating ρ(ψ, γ) over all specular angles ψ gives the diffuse intensity h(γ) in the direction γ. There are several natural requirements on ρ, including positivity and finite support, i.e., it should have nonzero values on a finite domain. In Figure 4b, its support may be considered the shading, where darker values represent a higher density, and the support is clearly a function of the angles, i.e., the scattering is spatially varying. Notice that the second energy balance in equation (7) is trivially fulfilled by direct substitution:(9a) (9b)which are the same after a change of integration order. Before defining this density ρ explicitly, we note that it is closely related to the bidirectional reflectance distribution function (BRDF), which is defined as the ratio between the outgoing radiance and the incoming irradiance on some small area element [15]. In a future publication, we shall explore this connection, and show for which circumstances the choice of ρ we make below is valid.
For now, we shall make an ad hoc choice of the density ρ. In particular, consider(10)where p is a function describing the redistribution of light, subject to an energy constraint we shall formulate momentarily. Physically, this choice can be motivated as follows. In the specular case, i.e., Figure 4a, p(a(ψ, γ); ψ) = p(γ − ψ) = δ(γ − ψ), where δ represents the Dirac delta function, meaning the light will be scattered in exactly one direction γ ≡ ψ. When p is some other appropriate function, light in direction ψ is scattered over multiple angles and we have a situation similar to that in Figure 4b, i.e., the choice of ρ in equation (10) represents the physical properties of light scattering. Note, however, that p includes the parameter ψ, which highlights the possibility of unique p functions for each specular ray, which is what is schematically shown in Figure 4b since the support varies with . Inserting this density in equation (8a), transforming the integration variable γ to α, and noting that with our maps, the Jacobian ∂s/∂α = 1, we get(11)where α_{1} = min{a(ψ, γ)ψ ∈ [ψ_{1}, ψ_{2}], γ ∈ [γ_{1}, γ_{2}]} and α_{2} = max{a(ψ, γ)ψ ∈ [ψ_{1}, ψ_{2}], γ ∈ [γ_{1}, γ_{2}]}. From equation (12) and using the fact that g(ψ) ≥ 0, it is clear that p(α; ψ) ≥ 0, i.e., with this choice of ρ, p becomes a probability density function (PDF).
2.1.3 Integral equation
Let us now focus on the second integral relation – equation (8b). Substituting our choice of ρ(ψ, γ) from equation (10) yields(12)
We once again utilise α = a(ψ, γ) to change the integration variable from ψ to α, together with ψ = s^{−1} (γ; α), to get(13)where α_{1} and α_{2} were defined after equation (13). Here, we note that this is a Fredholm integral equation of the first kind for g, as p depends on both γ (via ψ = s ^{−1} (γ; α)) and α. That is, p is a spatially varying kernel function, h is the prescribed target and g is to be determined.
Under the assumption of isotropic scattering, the explicit ψdependence in p is omitted, meaning we get p(a(ψ, γ)), or simply p(α). The ψ vs. γ plot for such a situation is shown in Figure 5. In contrast to Figure 4b, the support of ρ is now a band of constant width, and the data represent that of Example # in Section 4. Inserting a(ψ, γ) and s ^{−1}(γ; α) from equation (6) into equations (12) and (13) yields(14a) (14b)which are convolution integrals. We shall use the common notation of h(γ) = (p*g)(γ) for the convolution in equation (14a). Due to the commutativity property of convolution integrals, an equivalent definition is h(γ) = (g*p)(γ) in equation (14b) ([14], p. 309). Obtaining g is now a matter of deconvolving equation (14a) or (14b). There are a large number of different deconvolution methods, but not all are equally suitable for our purposes. In particular, g must be nonnegative within a finite domain. This can be achieved using an iterative ratio method, such as Gold’s method, or the more common Richardson–Lucy method of deconvolution [17]. We shall return to this topic in Section 4.
2.2 Rotationally symmetric reflectors
To introduce rotationally symmetric, consider Figure 6a. The situation is as follows. A pointsource is located at the origin . The reflector is parametrised by , where u(φ) > 0 is at least twice continuously differentiable and is the radial unit vector in spherical coordinates. The angle φ ∈ [0, π) is measured from the positive zaxis and the angle ϑ ∈ (−π, π) is measured from the positive xaxis. Consider now a ray in direction emitted from the point source at . Following its trajectory, it strikes the reflector at a point , where the unit normal to the surface (not shown) is given by . Like in the cylindrically symmetric case, we adopt the convention , i.e., the normal points towards the light source. From the LoR, equation (1), we get an expression for the reflected ray .
Figure 6 Rotationally symmetric reflectors exhibiting arbitrary (a), and inplane (b) scattering; point source at . 
2.2.1 Inplane scattering
If, in addition to the reflector being rotationally symmetric, the scattered direction is assumed to be in the plane of incidence, the problem may be analysed in two dimensions. The fact that the specular problem reduces to two dimensions is shown in [10] of Chapter 3, and inplane scattering preserves this symmetry. The situation is shown in Figure 6b, depicting the crosssection in the plane of incidence (here taken as the xzplane to distinguish from the xyplane used in the cylindrically symmetric problems). The twodimensional reflector is parametrised by , where , measured from the positive zaxis. The source ray intersects the reflector at a point with unit normal , and the specular direction , (not shown) measured from the positive zaxis, is given by the LoR, equation (1). The scattered direction , γ ∈ [0, π) measured from the positive zaxis, is given by a rotation by a stochastic parameter α of , in accordance with equation (2), around an axis parallel to the suppressed yaxis passing through .
Before proceeding, let us briefly discuss when such an approximation is valid. The rotationally symmetric reflector is selfexplanatory, but the inplane scattering is less straightforward. The proposed situation where this may hold is as follows. A rotationally symmetric reflector was machined in a manner which left the surface chiselled so that the chiselmarks follow a tightlywound spiral along the reflector. In this situation, we postulate that we may study the crosssection in the plane of incidence, since the symmetry of the problem is preserved. Within these restrictions, i.e., a rotationally symmetric reflector subject to our postulate of inplane scattering only, we may readily use the formulae in the prior sections, with a few changes that shall be highlighted shortly. We furthermore assume that the intensity distributions are rotationally symmetric, such that the following description will suffice [lm/rad]:

source intensity distribution f(φ) > 0, φ ∈ [φ_{1}, φ_{2}],

intermediate specular intensity distribution ,

diffuse target intensity distribution .
The energy balances in equation (7) become(15)where the sine terms come from the spherical area elements. Following the procedure in the previous section, let us introduce the density ρ(ψ, γ) ≥ 0, ψ ∈ [ψ_{1}, ψ_{2}], γ ∈ [γ_{1}, γ_{2}], such that(16a) (16b)
Comparing these equations to equation (8), notice that they are the same up to the sine terms from the spherical area elements. As such, the natural choice for ρ becomes (recall Eq. (12); assuming isotropic scattering, i.e., no explicit ψdependence)(17)
We now substitute this ρ into equation (16a) to get the normalisation condition(18)
Transforming γ to α yields(19)where α_{1} and α_{2} were defined after equation (11). The Jacobian ∂s/∂α = 1 via the mappings in equation (6). Whence, the normalisation of p is(20)
Substituting ρ, defined in equation (17), into equation (16b) yields(21)
Absorbing the sine terms in the intensity distributions by defining(22)and inserting a(ψ, γ) from equation (6), yields(23a) (23b)where the second equation is obtained by transforming ψ to α. Comparing these to equations (14a) and (14b), it is clear that they are the same, up to the sine terms from the spherical area elements in the modified distributions. As such, deconvolution is still a vital tool to obtain the specular target distribution used in the reflector design procedure.
3 Specular reflector design
Having obtained a specular target distribution via deconvolution, our goal is as follows. Determine a specular reflector which transforms the given source distribution f into the given target distribution g. The approach we have chosen involves solving two ordinary differential equations (ODEs) for the radius function u(φ) and the mapping m(φ), which together fully characterise the reflector. This is similar to the approach outlined in [13] of Chapter 3.3. Note that the resulting equations are scaleinvariant, so the physical size of the final reflector is arbitrary and can be decided by the optical engineer by choosing appropriate units of length.
3.1 Cylindrically symmetric reflectors
Recall that the reflector is parametrised by , and that φ is measured counterclockwise from the positive xaxis (refer to Fig. 2). To start, note that a tangent vector to the reflector is(24)where and are the standard unit vectors in polar coordinates, and where we made use of the relation . The corresponding normal vector can be constructed by rotating counterclockwise, i.e.,(25)where the rotation matrix R was initially defined in equation (2). The associated unit normal is(26)where we made use of the fact that and . Let , so that(27)
Let us compute(28)which we shall use momentarily. In doing so, we made use of the relation . Note that , indicating we rotated τ the correct way to get n. From the LoR, equation (1), we get(29)
By geometrical arguments – see Figure 2a – it is clear that , such that together with equation (28), we get(30)or, equivalently(31)where we used the tangent halfangle relation ([16], p. 127) and switched from ψ to m(φ) in the last step to highlight that this is indeed the specular map – recall Figure 3, which has an explicit φdependence. To solve this ODE, we shall make use of the arbitrary initial condition v(φ_{1}) = 0. We thus have the initial value problem (IVP)(32)
This IVP describes the rate of change of the reflector radius function u, via v = ln(u), such that each incident angle φ is converted into the corresponding specular angle ψ = m(φ). At this stage, m is still unknown, and it is determined by the first energy balance in equation (7).
As an example, we consider a monotonically increasing or decreasing optical map m(φ). Suppose we have a monotonically increasing function , together with the condition , such that the reflected rays do not intersect, i.e., the ray bundle is divergent. Since m_{div} is by construction a valid solution, i.e., it achieves g, given f, the following must hold for all φ ∈ [φ_{1}, φ_{2}] (recall that ψ_{1} < ψ < ψ_{2})(33)i.e., the total flux in the interval [φ_{1}, φ] emitted by the source is equivalent to the reflected flux in the image interval [ψ_{1}, m_{div}(φ)]. Differentiation with respect to φ immediately yields the IVP(34)
Analogous considerations with a monotonically decreasing function , where we instead have a convergent ray bundle and where , yield, for all ,(35)or in terms of the equivalent IVP,(36)
Once the desired m(φ) has been obtained by solving the corresponding IVP, it is substituted into equation (32), which is then solved, yielding and consequently , which fully characterises the reflector. Throughout this paper, we have utilised Matlab’s ode15s routine to solve the IPVs.
3.2 Rotationally symmetric reflectors
In the case of rotationally symmetric reflectors, we measure φ ∈ [0, π) from the positive zaxis, in the plane of incidence – recall Figure 6b. Since ϑ is constant, the reflector is parametrised by , where is the radial unit vector in this particular polar coordinate system. Thus, the tangent vector becomes(37)where is the angular unit vector in this polar coordinate system. To obtain the unit normal pointing towards the source, we rotate the tangent vector clockwise, i.e.,(38)
Following the procedure from the previous section, we end up with a unit normal(39)after introducing . Finally, we consider the geometry of the situation (see Fig. 6b; ψ is analogous to γ) to conclude that , so that together with the law of reflection, equation (1), the boundary condition v(φ_{1}) = 0, and the tangent halfangle relation, we once again arrive at the IVP in equation (32).
As in the cylindrically symmetric case, suppose we have a monotonically increasing specular map . Then, for any ,(40)or, formulated as an IVP(41)
Similarly, for a monotonically decreasing ,(42)and the IVP becomes(43)
With these changes in mind, and the knowledge that the IVP for v remains the same, we can conclude that the design procedure outlined previously may be used without further modifications.
3.3 Raytracing
Irrespective of how the reflector is computed, and whether it is cylindrically or rotationally symmetric, a validation method is required. The natural choice is raytracing, and we have written our own twodimensional raytracer, since we must have full control of how the scattering occurs at the surface in order to validate our model. In particular, source, specular and diffuse (i.e., scattered) rays are all collected. The source rays are generated from the appropriate source distribution using Matlab’s rand routine, followed by an intersection computation. When computing the reflector, we are left with discrete data points, and these form reflector segments via piecewiselinear interpolation. As such, all rays that fall within a reflector segment will use the same normal in the computation of the reflected direction. The intersection point for a given source ray is decided by computing the angle of a line connecting the centre of each reflector segment and the origin , and then comparing these angles to that of the generated ray using Matlab’s dsearchn routine. The specular rays are then computed using the law of reflection, equation (1), whilst the diffuse rays are computed using a rotation matrix – recall equation (2). The stochastic parameter α is sampled from the chosen scattering function p. We shall use either a Gaussian and Matlab’s randn routine or a Lorentzian (Cauchy distribution) and Matlab’s rand routine in the corresponding cumulative distribution function. The ray collection is performed by equidistantly dividing the relevant angular domain ((−π, π) or [0, π), for cylindrically and rotationally symmetric distributions, respectively), thus forming collection bins. The centres of the collection bins are known, and the appropriate bin for a given ray is then computed via a nearest point search using dsearchn, and the number of rays collected by that bin is incremented by unity. After this, the process is repeated until the requested number of rays have been traced though the system. Finally, the number of rays per collection bin is converted into an intensity by dividing the probability of falling in each bin by the size of the bins and multiplying with the total flux of the source, i.e.,(44)for cylindrically symmetric reflectors, and for the jth bin. Here, 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 Δφ is the angular size of the collection bins. The total flux of the source is given by the integral over f. For rotationally symmetric reflectors, we have(45)
More details regarding raytracing are available in [18], p. 34.
Angle convention
We adopt the (−π, π) angle convention for cylindrically symmetric reflectors, so that we can make use of Matlab’s atan2 function to compute the angles of the rays. The rotationally symmetric problems will still use the polar [0, π) convention. In addition to this, we shall define our target distributions in terms of κ(ψ) and κ(γ), where(46)
Validation criterion
We utilise raytracing as a validation technique, so let us define a criterion to quantify the differences between the raytraced and predicted distributions. Specifically, we use the root mean square (RMS) error defined as follows(47)where N is the number of collection bins of h*, and the star indicates raytraced distributions. In most of our examples, we use the deconvolved (subscript “dc”) specular distribution g_{dc}, obtained by deconvolving equations (14) or (23) when designing the reflectors. In this case, h and h* will be h_{rc} and , respectively, where the “rc” subscript signifies “reconvolution”, that is .
4 Examples
This section presents three sample problems: two exhibiting cylindrical symmetry and one with rotational symmetry. To verify our model, we prescribe the specular target distribution g exactly, and construct the diffuse target distribution h by convolving g and the chosen scattering PDF p. We then compute the deconvolved specular distribution g_{dc}, which is used to design the reflectors. Finally, we raytrace the system and compare the result to our prediction. In the rotationally symmetric example, we no longer know the exact g, but rather we prescribe h exactly. This is more similar to how we envision an optical designer working with our model.
The computation time for these examples may be divided as follows (times given for a portable workstation with 4 cores, 8 threads @ 4 GHz and 32 GB of RAM):

Computing the deconvolution:

Depends on the number of distribution bins.

Can be done very efficiently using Matlab routines, such as the builtin deconvolucy routine we used.

Using efficient routines, this step takes a few seconds on our machine.


Computing the specular reflector:

Depends on the number of reflector segments, distribution bins, and the stiffness of the ODEs.

Involves solving two IVPs, which can be done using various methods, such as Matlab’s ode15s routine, which we used.

For simple problems, like our examples, this step takes up to a minute; for difficult problems it can take a few minutes.


Validating the reflectors using raytracing:

Depends on the number of rays traced.

This step can be parallelised to a very large extent, and modern consumer GPUs have dedicated raytracing hardware.

Using our raytracer in Matlab without any raytracing hardware and tracing a total of rays, this step takes a few minutes.

From the above list, it is clear that the step we propose in this paper, the deconvolution step, takes the least amount of time.
4.1 Example #1: smooth target distribution
The specular problem consists of a homogeneous source f being transformed into two partly overlapping Gaussians g. As for the choice of p, we opted for a Gaussian centred around α = 0°, with standard deviation σ = 10°. This is supposed to represent relatively minor scattering when compared to, e.g., Lambert’s cosine law, whilst still being a significant deviation from a specular reflector. The diffuse distribution h is the convolution between p and g, i.e., h = p*g. Worth noting is that Gaussians do not have finite support, meaning we need to truncate the nonzero values outside of [α_{1}, α_{2}] when performing the (de)convolution. We renormalised p after truncation to ensure that . The problem is summarised in the box below, where(48)represents the Gaussian, centred at μ with standard deviation σ. The value of φ_{2} was chosen such that energy is conserved up to .
Example #1: Smooth target distribution
φrange: [φ_{1}, φ_{2}] = [−π/4, 29π/75]
ψrange: [ψ_{1}, ψ_{2}]: see text
αrange: [α_{1}, α_{2}] = [−π, π]
Source distribution:
Specular target distribution:
Surface scattering function:
Diffuse distribution prediction:
Boundary condition: u(φ_{1}) = 1
Since we are interested in validating the whole proposed solution method, the first step is to find the specular target distribution g_{dc} by deconvolution. We shall utilise Richardson–Lucy deconvolution, and in particular Matlab’s deconvlucy routine. This iterative ratio method has numerous benefits compared to direct methods, most crucial for our purposes being guaranteed positivity of the solution. The functions f, g, g_{dc}, p, h, and h_{rc} are shown in Figure 7, for deconvlucy’s default settings of 10 iterations with 128 sampling points. Clearly, the recovered g_{dc} resembles the original g relatively well (and it could be improved by increasing the number of deconvolution iterations), and the reconvolved h_{rc} = g_{dc}*p is nearly identical to h.
Figure 7 Distributions in Example #1; 128 sample points. 
The next step is to design the reflectors. In order to use the procedure outlined in Section 3, we need the limits ψ_{1} and ψ_{2}. Recall that these should represent the support of g (or, rather, g_{dc} in this case). In this example, the limits are ambiguous due to the Gaussians. We computed the limits by fixing a threshold η = 0.001 and locating the two extrema of κ(ψ) where g_{dc} (ψ) = η, using piecewiselinear interpolation between the data points of g_{dc}. The limits are shown in Figure 8, and the reflectors are shown in Figure 9. In this case, we do not know the exact solutions, so we shall not discuss the reflectors further for this example. We note that one could renormalise g_{dc} to ensure more accurate energy conservation, but this has not been done in the data shown. We have used 1024 sample points for the reflectors in an attempt to minimise discretisation errors due to the reflectors when validating the m_{conv} reflector using raytracing (recall that each ray striking a reflector segment uses the same normal). The sample points are equidistant in the [φ_{1}, φ_{2}]range. The raytraced distributions and the RMS error from equation (47) are shown in Figure 10, where we see that the source sampling is appropriate, and the resulting distributions are well predicted by our model. In addition, the convergence shows the expected behaviour of Monte Carlo raytracing ([18], p. 9), where N_{ r } is the number of rays traced through the system.
Figure 8 The boundaries used as the support of g_{dc} in Example #1. 
Figure 10 Example #1; raytraced distributions of m_{conv} reflector with g_{dc}; 10^{6} rays. 
There are a couple of regions where the raytraced distributions deviate from our predictions. Specifically, near and , near and near both peaks of the Gaussians. The discrepancies in are due to the bins not aligning perfectly with the support of , such that part of a collection bin may cross the φ_{1} and φ_{2}boundaries. The discrepancies of are presumably partly due to energy not being perfectly conserved, and partly from the discretisation of the reflector surface, in addition to the aforementioned binning issue. Additionally, the very peaks of the Gaussians are only one or two data points wide, and achieving that level of precision is no easy feat, using a numerical scheme. Keeping all of these factors in mind, the results are promising, and from the RMS error plot, we see that increasing the number of rays is likely to improve the result further.
4.2 Example #2: Block function as target distribution
We now move on to our second example, which at first appears much simpler, but will prove to be quite a challenge for our numerical scheme. The specular problem consists of homogeneous illumination of a circular disk within and is homogeneous on . The scattering PDF is still a Gaussian, this time with a standard deviation of σ = 5°. We employ a similar approach to the first example, i.e., prescribe g, compute and attempt to recover g via deconvolution, then validate the reflectors we designed using with raytracing. The example is outlined in the box below. We wish to highlight that the density for this example is shown in Figure 5.
Example #2: Block function as target distribution
φrange: [φ_{1}, φ_{2}] = [−π/4, π/4]
ψrange: [κ(ψ_{1}), κ(ψ_{2})] = [−π/4, π/4]
αrange: [α_{1}, α_{2}] = [−π, π]
Source distribution:
Specular target distribution:
Surface scattering function:
Diffuse distribution prediction:
Boundary condition: u(φ_{1}) = 1
Using the default settings of deconvlucy (10 iterations) yields g_{dc} in Figure 11, where we immediately see that it deviates significantly from the original g. Readers who are familiar with signal theory are likely not surprised by this, as representing a block function in Fourier space requires an infinite number of frequencies. Let us attempt to increase the number of deconvolution iterations by an order of magnitude – see Figure 12a. This shows a slight improvement, but we are still quite far from the original g. As such, let us further increase the number of iterations by two orders of magnitude to get the result in Figure 12b, which is certainly a lot closer to the original g. The RMS error ε(g, g_{dc}), recall equation (47), decreased from 0.055 to 0.020 and 0.005, for 10, 10^{2} and 10^{4} iterations, respectively. In a real problem, this metric would not be available to us, so we would have to rely on the RMS error ε(h, h_{ rc }), which decreased from 10^{−3} to 10^{−4} and 10^{−6}. Based solely on ε(h, h_{rc}), it is not unreasonable that one might design the reflector using the first g_{dc}, so we shall include it as a worstcase scenario, as well as the best g_{dc,}, in the sense that it has the lowest RMS error.
Figure 11 Initial distributions in Example #2; 128 sample points. 
Figure 12 Example #2; deconvlucy deconvolution with 10^{2} iterations (a), and 10^{4} iterations (b). 
Turning to the topic of reflector design, consider the reflectors in Figures 13a and 13b, designed using the original g and the deconvolved g_{dc} in Figure 11, respectively. We note that the exact solutions to this problem are given in [13] (p. 28) as a circle segment and a straight line, i.e., we recover them using our numerical scheme. To the naked eye, the two figures appear nigh identical, and it is only when we plot the difference in radii of the m_{conv} reflectors in Figure 14 that we can appreciate the differences. From the raw (or unaltered) graph, we postulate that the deviations can be decomposed into a sloped straight line and comparatively small oscillations. In order to better appreciate the oscillations, we thus subtracted a linear correction factor from the raw data. This reveals the profile of g_{dc} (ψ), ψ ∈ [ψ_{1}, ψ_{2}], present in the reflector surface. Notice that the order of magnitude of the oscillations in the reflector surface is significantly smaller than that of the oscillations in the specular target distribution. In other words, small changes in the reflector can result in significant changes to the resulting distribution.
Figure 13 Example #2; reflectors designed using g in Figure 11 (a), and g_{dc} after 10 deconvolution iterations in Figure 11 (b); 1024 sample points. 
Figure 14 Difference in reflector radii of the m_{conv} reflectors in Figures 13a and 13b; slope of the subtracted linear correction term in the right panel was 4.31 × 10^{−6}. 
Shifting focus to raytracing, the results are shown in Figures 15a–15c with g and g_{dc} from Figure 11 (10 iterations) and g_{dc} from Figure 12 (10^{4} iterations), respectively. In all these cases, we used a total of 10^{6} rays, and the m_{conv} reflectors. It is clear that the scheme works well, and the issues we see were explained when discussing the previous example. That is, binning and numerical errors due to discretisation. A slight asymmetry appears in the results for the reflectors designed using the deconvolved g_{dc} distributions. This can also be seen from the slope in Figure 14, so the cause appears to be somewhere in the numerical computation of the reflectors, presumably due to integration from left to right. The discrepancy is minor, so an attempt to correct it has not been made, as it is clear that the model predicts the scattered distribution very well, and that the resulting reflector can be validated using raytracing.
Figure 15 Example #2; raytraced distributions of m_{conv} reflectors with g from Figure 11 (a), g_{dc} from Figure 17 (b), and g_{dc} from Figure 12b (c); 10^{6} rays. 
4.3 Example #3: Lorentzian scattering function
The rotationally symmetric example we have chosen differs from the previous examples in two major ways. The first is that we no longer know the exact g, but rather we prescribe an exact h. This is more similar to how we envision an optical designer would incorporate our model in their workflow. The second difference is that the scattering function is a Lorentzian (also known as a Cauchy distribution). This is significant for two reasons, namely that machined mirrors often exhibit this type of bidirectional reflectance distribution function (BRDF) ([19], Chap. 4), and because the tails fall to zero at a significantly lower rate. The latter fact means that more largeangle scattering will occur when compared to the Gaussians we have used thus far. As such, we increase the relevance of the method whilst testing the limits of our model. The example is outlined in the box below, where(49)is the Lorentzian function with a full width at half maximum (FWHM) of 2σ; σ is often denoted γ in literature, but not here for obvious reasons. Note that we again truncated the values of p outside of [α_{1}, α_{2}] and renormalised, such that .
Example #3: Lorentzian Scattering Function
φrange: [φ_{1}, φ_{2}] = [π/4, 3π/4 − 0.34]
γrange: [κ(γ_{1}), κ(γ_{2})] = [1.015, 2.222]
αrange: [α_{1}, α_{2}] = [−π, π]
Source distribution:
Diffuse target distribution:
Surface scattering function:
Boundary condition: u(φ_{1}) = 1
The distributions are shown in Figure 16, where we have opted to absorb the sine terms from the energy balances, equation (15), into the distributions (indicated by the tilde). The deconvolved specular target distribution was computed using the default deconvlucy settings with 10 iterations. Before designing the reflectors, let us briefly consider what to expect from the final raytraced distributions. In particular, since we have prescribed h, rather than g, we are no longer guaranteed that the deconvolution converges. We can get an appreciation for this by comparing the reconvolved with our prescribed in Figure 16. This may seem like a disappointing result, but let us compare , and , where the latter would be the diffuse result if we disregarded scattering in the design procedure entirely, i.e., if we designed the reflectors using f and took as g in the design procedure, and raytraced the optical system using our scattering model. All of these distributions are shown in Figure 17, and the RMS error , whilst , i.e., our approach represents an improvement of approximately 40%. Visually, we see that the problematic regions for are partly the peaks and partly near and . The deviation close to the peaks could perhaps be improved by increasing the number of deconvolution iterations, but the problems close to the boundaries are not solvable in our model. This is due to an inherent “maximum steepness” dictated by the least steep function we are deconvolving (in this case p), and it is a property of deconvolution. The only way to reach the correct values at the boundaries is to introduce negative values outside of them, which would be unphysical.
Figure 16 Initial distributions in Example #3; 128 sample points. 
Figure 17 The prescribed and predicted diffuse targets in Example #3. 
We are now ready to design the reflectors. By fixing η = 0.001 and locating the two values of where , we found and , see Figure 18. The reflectors computed using g_{dc} as the specular target distribution and f as the source distribution are shown in Figures 19a and 19b, where the latter is a threedimensional version of the m_{conv} reflector. The raytraced distributions are shown in Figure 20a, where we see that the source sampling is correct, as is the resulting diffuse distribution . As for the intermediate specular target distribution, some deviations from the target may be seen, especially near the first peak. These deviations are presumably due to difficulties solving the relevant IVPs that give the reflector radius function u, and they are likely the reason why our RMS error convergence slows down after approximately rays. For the sake of completeness, we also raytraced the reflector in Figure 20b. Here, we only show , , and , since the source sampling is identical. It is clear that the diffuse distribution is still very closely achieved, whilst the specular distribution deviates in more obvious ways than with the m_{conv} reflector. This highlights the nontriviality of designing specular reflectors, rather than any apparent flaw in our model of scattering, and it could perhaps be improved by increasing the sampling frequency of the distributions, or altering the number of samples in the reflectors themselves. This has not been attempted, since we are mostly interested in validating our model for scattering, and indeed, even with these deviations from the specular target distribution, the effect of scattering smoothes them out greatly.
Figure 18 The κ(ψ)boundaries used as the support of g_{dc} in Example #3. 
Figure 19 Example #3; reflectors designed using g_{dc} in Figure 16 (a), and a threedimensional version of the m_{conv} reflector (b). 
Figure 20 Example #3; raytraced distributions of m_{conv} reflector with g_{dc} (a), and m_{div} reflector with g_{dc} (b); 10^{6} rays. 
5 Conclusions
We have presented a novel modelling approach to include surface scattering in the design of reflectors as part of optical systems. The approach is inspired by concepts from optimal mass transport theory, and it relies on energy conservation. In the case of isotropic inplane scattering and cylindrical or rotational symmetry, the forward expression of the scattered light reduces to a convolution integral between a probability density function and a specular target function. By prescribing a desired diffuse target distribution, and the scattering probability density function, one can solve for the specular target distribution using deconvolution methods from literature, and then compute the reflectors using purely specular design procedures. As such, including the effects of scattering can be considered a preprocessing step, and all the mature specular reflector design procedures remain essential. This gives the optical designer a greater ability to use scattering to their advantage, or mitigate it in applications where it is undesirable.
We view this work as a proof of concept and as an important first step towards a generalised threedimensional model of surface scattering using geometrical optics, with the ultimate goal of computing freeform reflectors with scattering surfaces. Lifting the requirement of inplane scattering would significantly increase the applicability of our model and we are actively working in this direction. Additionally, a manuscript exploring the connection between this work and BRDF functions is currently in review.
Funding
This work was partially supported by the Dutch Research Council (Dutch: Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)) through grant P1536.
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.
References
 Hecht E. (2017) Optics, 5th edn., Pearson Education, Incorporated. [Google Scholar]
 Strutt J.W. (1907) On the dynamical theory of gratings, Proc. R. Soc. Lond. A 79, 399–416. [NASA ADS] [CrossRef] [Google Scholar]
 Rice S.O. (1951) Reflection of electromagnetic waves from slightly rough surfaces, Commun. Pure Appl. Math. 4, 2–3, 351–378. [CrossRef] [Google Scholar]
 Beckmann P., Spizzichino A. (1987) The scattering of electromagnetic waves from rough surfaces, Artech House Radar Library, Artech House. [Google Scholar]
 Harvey J.E. (2019) Understanding surface scatter phenomena: a linear systems formulation, SPIE Press monograph, SPIE Press. [Google Scholar]
 Church E.L., Zavada J.M. (1975) Residual surface roughness of diamondturned optics, Appl. Opt. 14, 8, 1788–1795. [NASA ADS] [CrossRef] [Google Scholar]
 Harvey J.E. (1976) Lightscattering characteristics of optical surfaces, PhD thesis, University of Arizona. [Google Scholar]
 Krywonos A., Harvey J.E., Choi N. (2011) Linear systems formulation of scattering theory for rough surfaces with arbitrary incident and scattering angles, J. Opt. Soc. Am. A 28, 6, 1121–1138. [NASA ADS] [CrossRef] [Google Scholar]
 McNamara D.A., Pistorius C.W.I., Malherbe J.A.G. (1990) Introduction to the uniform geometrical theory of diffraction, Antennas and Propagation Library, Artech House. [Google Scholar]
 Villani C. (2003) Topics in optimal transportation, Graduate Studies in Mathematics. American Mathematical Society. [Google Scholar]
 Prins C.R. (2014) Inverse methods for illumination optics, PhD thesis, Eindhoven University of Technology. [Google Scholar]
 Romijn L.B., ten Thije Boonkkamp J.H.M., IJzerman W.L. (2020) Inverse reflector design for a point source and farfield target, J. Comput. Phys. 408, 109283. [NASA ADS] [CrossRef] [Google Scholar]
 Maes M.J.J.J.B. (1997) Mathematical methods for reflector design, PhD thesis, University of Amsterdam. [Google Scholar]
 Lin R.J., Tsai M.S., Sun C.C. (2015) Novel optical lens design with a light scattering freeform inner surface for LED down light illumination, Opt. Exp. 23, 13, 16715. [NASA ADS] [CrossRef] [Google Scholar]
 Nicodemus F.E. (1965) Directional reflectance and emissivity of an opaque surface, Appl. Opt. 4, 7, 767–775. [NASA ADS] [CrossRef] [Google Scholar]
 Råde L., Westergren B. (2004) Mathematics handbook for science and engineering, SpringerVerlag, Berlin Heidelberg. [CrossRef] [Google Scholar]
 Jansson P.A. (1997) Deconvolution of images and spectra, Academic Press. [Google Scholar]
 Filosa C. (2018) Phase space ray tracing for illumination optics, PhD thesis, Eindhoven University of Technology. [Google Scholar]
 Stover J.C. (2012) Optical scattering: measurement and analysis, Press Monograph Series, SPIE Press. [Google Scholar]
All Figures
Figure 1 The inverse problem of reflector design with a scattering surface; the first step (deconvolution) yields the specular target distribution, which in turn can be used to compute the reflector by solving the inverse specular problem. The data is taken from Example #1 in Section 4. 

In the text 
Figure 2 Reflectors exhibiting specular reflection (a), and light scattering (b); linesource along suppressed zaxis passing through . 

In the text 
Figure 3 Relations between unit vectors and angles; all symbols are defined in the text. 

In the text 
Figure 4 Schematic maps for specular reflectors (a), and diffuse reflectors (b). 

In the text 
Figure 5 Diffuse map (isotropic scattering; example # in Sect. 4). 

In the text 
Figure 6 Rotationally symmetric reflectors exhibiting arbitrary (a), and inplane (b) scattering; point source at . 

In the text 
Figure 7 Distributions in Example #1; 128 sample points. 

In the text 
Figure 8 The boundaries used as the support of g_{dc} in Example #1. 

In the text 
Figure 9 Example #1; reflectors designed using g_{dc} in Figure 7; 1024 sample points. 

In the text 
Figure 10 Example #1; raytraced distributions of m_{conv} reflector with g_{dc}; 10^{6} rays. 

In the text 
Figure 11 Initial distributions in Example #2; 128 sample points. 

In the text 
Figure 12 Example #2; deconvlucy deconvolution with 10^{2} iterations (a), and 10^{4} iterations (b). 

In the text 
Figure 13 Example #2; reflectors designed using g in Figure 11 (a), and g_{dc} after 10 deconvolution iterations in Figure 11 (b); 1024 sample points. 

In the text 
Figure 14 Difference in reflector radii of the m_{conv} reflectors in Figures 13a and 13b; slope of the subtracted linear correction term in the right panel was 4.31 × 10^{−6}. 

In the text 
Figure 15 Example #2; raytraced distributions of m_{conv} reflectors with g from Figure 11 (a), g_{dc} from Figure 17 (b), and g_{dc} from Figure 12b (c); 10^{6} rays. 

In the text 
Figure 16 Initial distributions in Example #3; 128 sample points. 

In the text 
Figure 17 The prescribed and predicted diffuse targets in Example #3. 

In the text 
Figure 18 The κ(ψ)boundaries used as the support of g_{dc} in Example #3. 

In the text 
Figure 19 Example #3; reflectors designed using g_{dc} in Figure 16 (a), and a threedimensional version of the m_{conv} reflector (b). 

In the text 
Figure 20 Example #3; raytraced distributions of m_{conv} reflector with g_{dc} (a), and m_{div} reflector with g_{dc} (b); 10^{6} rays. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.