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 so-called 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 trial-and-error 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 wide-angle 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 moderate-angle 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 Harvey-Shack 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 in-plane 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 pre-processing 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 two-dimensional 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.

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 pre-processing step, such as in the algorithm below.

The structure of this paper is as follows. Cylindrically symmetric problems with in-plane scattering are covered first, together with a few words about deconvolution in Section 2.1, followed by rotationally symmetric problems with in-plane scattering in Section 2.2. Next, two-dimensional 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.

Starting with cylindrically symmetric problems, consider the situation depicted in Figure 2. This specular problem can be viewed as a cross-section of a translationally invariant problem, such as an extruded optical element, with a line-source along the suppressed z-axis. For such a system, the specular problem may be analysed in two-dimensions ([10], Chap. 3). Thus, the intensities and reflector surfaces are independent of z and we may study a cross-section in the plane of incidence, which is spanned by the source and reflected rays, along ŝ and t̂, respectively, and which contains the unit normal n̂. The hat (̂) indicates unit vectors throughout this manuscript, and positive angles are by convention counter-clockwise. The reflector is parametrised by r(φ)=u(φ)êr, where u(φ) > 0 is at least twice continuously differentiable and êr=(cos(φ),sin(φ))ú is the radial unit vector in polar coordinates. The angle φ is measured counter-clockwise from the positive x-axis, and it fully characterises the source ray along ŝ≡êr, emitted form the line source at the origin O. The source ray along ŝ intersects the reflector at some point P, where the unit normal of the reflector is given by n̂. We take the convention ŝ⋅n̂<0, 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 t̂=(cos(ψ),sin(ψ))ú, i.e.,(1)t̂=ŝ-2(ŝ⋅n̂)n̂,where we denote the angle between the positive x-axis and t̂ by ψ.

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 cross-section of the translationally invariant problem. This assumption is made in order to preserve the two-dimensional 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 line-source parallel to the extrusion direction. This is because the process of extrusion leads to grooves along the extrusion direction, so that the cross-sectional 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 two-dimensional surface scattering model. Here, the source ray along ŝ gets mapped to a scattered ray along û=(cos(γ), sin(γ))T, where γ is measured counter-clockwise from the positive x-axis. The scattered direction û can be described as a rotation of t̂ by a stochastic parameter α around the axis parallel to the z-axis passing through P, i.e.,(2)û=R(α)t̂, where R(α):=cos(α)-sin(α)sin(α)cos(α).

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.

To express reflection and scattering in terms of angles, let us introduce two mappings which give the reflected and scattered directions, i.e.,(3)m(φ)=ψ and s(ψ;α)=γ,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 three-dimensional case, which we intend to treat in a future publication. In addition to these maps, we require their inverses to exist,(4)m-1(ψ)=φ and s-1(γ;α)=ψ.

Finally, we define a mapping yielding the scattering angle α, for fixed ψ and γ, i.e.,(5)a(ψ,γ)=α.

For a schematic summary, see Figure 3.

Explicitly, the maps are(6)m(φ)=φ+arccos(1-2(ŝ⋅n̂)2)≡ψ,m-1(ψ)=ψ-arccos(1-2(ŝ⋅n̂)2)≡φ,s(ψ;α)=ψ+α≡γ,s-1(γ;α)=γ-α≡ψ,a(ψ,γ)=γ-ψ≡α,where the first two relations follow from cos(ψ-φ)=ŝ⋅t̂ – 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.

Having presented the mappings for the angles, we introduce the intensity distributions (illuminance) [lm/rad]:

source intensity distribution f(φ) > 0, φ ∈ [φ₁, φ₂],

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 [Γ1,Γ2]⊂[γ1,γ2].

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, ψ ∈ [ψ₁, ψ₂], γ ∈ [γ₁, γ₂], with properties(8a)∫γ1γ2ρ(ψ,γ) dγ=g(ψ),(8b)∫ψ1ψ2ρ(ψ,γ) dψ=h(γ),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)∫γ1γ2h(γ) dγ=∫γ1γ2∫ψ1ψ2ρ(ψ,γ) dψdγ,(9b)∫ψ1ψ2g(ψ) dψ=∫ψ1ψ2∫γ1γ2ρ(ψ,γ) dγdψ,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)ρ(ψ,γ)=p(a(ψ,γ);ψ)g(ψ),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)∫α1α2p(α;ψ) dα=1,where α₁ = min{a(ψ, γ)|ψ ∈ [ψ₁, ψ₂], γ ∈ [γ₁, γ₂]} and α₂ = max{a(ψ, γ)|ψ ∈ [ψ₁, ψ₂], γ ∈ [γ₁, γ₂]}. 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).

Let us now focus on the second integral relation – equation (8b). Substituting our choice of ρ(ψ, γ) from equation (10) yields(12)h(γ)=∫ψ1ψ2p(a(ψ,γ);ψ)g(ψ) dψ.

We once again utilise α = a(ψ, γ) to change the integration variable from ψ to α, together with ψ = s⁻¹ (γ; α), to get(13)h(γ)=∫α1α2p(α;s-1(γ;α))g(s-1(γ;α)) ∂s-1(γ;α)∂α dα,where α₁ and α₂ 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⁻¹ (γ; α)) 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 #2 in Section 4. Inserting a(ψ, γ) and s⁻¹(γ; α) from equation (6) into equations (12) and (13) yields(14a)h(γ)=∫ψ1ψ2p(γ-ψ)g(ψ) dψ,(14b)h(γ)=∫α1α2p(α)g(γ-α) dα,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.

To introduce rotationally symmetric, consider Figure 6a. The situation is as follows. A point-source is located at the origin O. The reflector is parametrised by r(φ,ϑ)=u(φ)êr, where u(φ) > 0 is at least twice continuously differentiable and êr=(sin(φ)cos(ϑ),sin(φ)sin(ϑ),cos(φ))ú is the radial unit vector in spherical coordinates. The angle φ ∈ [0, π) is measured from the positive z-axis and the angle ϑ ∈ (−π, π) is measured from the positive x-axis. Consider now a ray in direction ŝ≡êr emitted from the point source at O. Following its trajectory, it strikes the reflector at a point P, where the unit normal to the surface (not shown) is given by n̂. Like in the cylindrically symmetric case, we adopt the convention ŝ⋅n̂<0, i.e., the normal points towards the light source. From the LoR, equation (1), we get an expression for the reflected ray t̂.

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 in-plane scattering preserves this symmetry. The situation is shown in Figure 6b, depicting the cross-section in the plane of incidence (here taken as the xz-plane to distinguish from the xy-plane used in the cylindrically symmetric problems). The two-dimensional reflector is parametrised by r(φ)=u(φ)êr, where êr=(sin(φ),cos(φ))ú, φ∈[0,π) measured from the positive z-axis. The source ray ŝ≡êr intersects the reflector at a point P with unit normal n̂, and the specular direction t̂=(sin(ψ),cos(ψ))ú, ψ∈[0,π) (not shown) measured from the positive z-axis, is given by the LoR, equation (1). The scattered direction û=(sin(γ),cos(γ))ú, γ ∈ [0, π) measured from the positive z-axis, is given by a rotation by a stochastic parameter α of t̂, in accordance with equation (2), around an axis parallel to the suppressed y-axis passing through P.

Before proceeding, let us briefly discuss when such an approximation is valid. The rotationally symmetric reflector is self-explanatory, but the in-plane scattering is less straight-forward. 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 chisel-marks follow a tightly-wound spiral along the reflector. In this situation, we postulate that we may study the cross-section 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 in-plane 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, φ ∈ [φ₁, φ₂],

The energy balances in equation (7) become(15)∫φ1φ2f(φ)sin(φ) φ=∫ψ1ψ2g(ψ)sin(ψ) ψ=∫γ1γ2h(γ)sin(γ) γ,where the sine terms come from the spherical area elements. Following the procedure in the previous section, let us introduce the density ρ(ψ, γ) ≥ 0, ψ ∈ [ψ₁, ψ₂], γ ∈ [γ₁, γ₂], such that(16a)∫γ1γ2ρ(ψ,γ)sin(γ) dγ=g(ψ),(16b)∫ψ1ψ2ρ(ψ,γ)sin(ψ) dψ=h(γ).

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)ρ(ψ,γ)sin(γ)=p(a(ψ,γ))g(ψ).

We now substitute this ρ into equation (16a) to get the normalisation condition(18)∫γ1γ2p(a(ψ,γ)) dγ=1.

Transforming γ to α yields(19)∫α1α2pα∂s(ψ; α)∂αdα=1,where α₁ and α₂ were defined after equation (11). The Jacobian |∂s/∂α| = 1 via the mappings in equation (6). Whence, the normalisation of p is(20)∫α1α2p(α)dα=1.

Substituting ρ, defined in equation (17), into equation (16b) yields(21)h(γ)sin(γ)=∫ψ1ψ2p(a(ψ,γ))g(ψ)sin(ψ)d ψ.

Absorbing the sine terms in the intensity distributions by defining(22)h̃(γ):=h(γ)sin(γ), g̃(ψ):=g(ψ)sin(ψ),and inserting a(ψ, γ) from equation (6), yields(23a)h̃(γ)=∫ψ1ψ2p(γ-ψ)g̃(ψ) dψ,(23b)h̃(γ)=∫α1α2p(α)g̃(γ-α)dα,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.

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 scale-invariant, so the physical size of the final reflector is arbitrary and can be decided by the optical engineer by choosing appropriate units of length.

Recall that the reflector is parametrised by r(φ)=u(φ)êr, and that φ is measured counter-clockwise from the positive x-axis (refer to Fig. 2). To start, note that a tangent vector to the reflector is(24)τ=r′(φ)=u′(φ)êr+u(φ)êφ,where êr=(cos(φ),sin(φ))ú and êφ=(-sin(φ),cos(φ))ú are the standard unit vectors in polar coordinates, and where we made use of the relation êr′=êφ. The corresponding normal vector can be constructed by rotating τ counter-clockwise, i.e.,(25)n=R(π/2)τ, R(π/2)=0-110,where the rotation matrix R was initially defined in equation (2). The associated unit normal is(26)n̂=n|n|=n|τ|=-u(φ)êr+u′(φ)êφu(φ)2+u′(φ)2,where we made use of the fact that R(π/2)êr=êφ and R(π/2)êφ=-êr. Let v(φ):=ln(u(φ)), so that(27)n̂=-êr+v′(φ)êφ1+v′(φ)2.

Let us compute(28)ŝ⋅n̂=-11+v′(φ)2,which we shall use momentarily. In doing so, we made use of the relation ŝ≡êr. Note that ŝ⋅n̂<0, indicating we rotated τ the correct way to get n. From the LoR, equation (1), we get(29)ŝ⋅t̂=1-2(ŝ⋅n̂)2.

By geometrical arguments – see Figure 2a – it is clear that ŝ⋅t̂=-ŝ⋅-t̂=cos(ψ-φ), such that together with equation (28), we get(30)cos(ψ-φ)=1-21+v′(φ)2,or, equivalently(31)v′(φ)=±1+cos(ψ-φ)1-cos(ψ-φ)=cotm(φ)-φ2,where we used the tangent half-angle 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(φ₁) = 0. We thus have the initial value problem (IVP)(32)v′φ=cotmφ-φ2v(φ1)=0, φ1<φ<φ2.

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 m(φ)=:mdiv(φ), together with the condition mdiv(φ1)=ψ1, 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 φ ∈ [φ₁, φ₂] (recall that ψ₁ < ψ < ψ₂)(33)∫φ1φf(φ̃) dφ̃=∫ψ1mdiv(φ)g(ψ̃) dψ̃,i.e., the total flux in the interval [φ₁, φ] emitted by the source is equivalent to the reflected flux in the image interval [ψ₁, m_div(φ)]. Differentiation with respect to φ immediately yields the IVP(34)m'div(φ)=f(φ)g(mdiv(φ)),mdiv(φ1)=ψ1.

Analogous considerations with a monotonically decreasing function m(φ):=mconv(φ), where we instead have a convergent ray bundle and where mconv(φ1)=ψ2, yield, for all φ∈[φ1,φ2],(35)∫φ1φf(φ̃) dφ̃=∫mconv(φ)ψ2g(ψ̃) dψ̃,or in terms of the equivalent IVP,(36)m′conv(φ)=-f(φ)g(mconv(φ))mconv(φ1)=ψ2, φ1<φ<φ2.

Once the desired m(φ) has been obtained by solving the corresponding IVP, it is substituted into equation (32), which is then solved, yielding v(φ) and consequently u(φ)=ev(φ), which fully characterises the reflector. Throughout this paper, we have utilised Matlab’s ode15s routine to solve the IPVs.

In the case of rotationally symmetric reflectors, we measure φ ∈ [0, π) from the positive z-axis, in the plane of incidence – recall Figure 6b. Since ϑ is constant, the reflector is parametrised by r(φ)=u(φ)êr, where êr=(sin(φ),cos(φ))ú is the radial unit vector in this particular polar coordinate system. Thus, the tangent vector becomes(37)τ=r′(φ)=u′(φ)êr+u(φ)êφ,where êφ=(cos(φ),-sin(φ))ú 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)n=R(-π/2)τ, R(-π/2)=01-10.

Following the procedure from the previous section, we end up with a unit normal(39)n̂=-êr+v′(φ)êφ1+v′(φ)2,after introducing v(φ):=ln(u(φ)). Finally, we consider the geometry of the situation (see Fig. 6b; ψ is analogous to γ) to conclude that ŝ⋅t̂=-ŝ⋅-t̂=cos(ψ-φ), so that together with the law of reflection, equation (1), the boundary condition v(φ₁) = 0, and the tangent half-angle 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 m(φ)=mdiv(φ). Then, for any φ∈[φ1,φ2],(40)∫φ1φf(φ̃)sin(φ̃) dφ̃=∫ψ1mdiv(φ)g(ψ̃)sin(ψ̃) dψ̃,or, formulated as an IVP(41)mdiv′(φ)=f(φ)sin(φ)g(mdiv(φ))sin(mdiv(φ))mdiv(φ1)=ψ1., φ1<φ<φ2.

Similarly, for a monotonically decreasing m(φ)=mconv(φ),(42)∫φ1φf(φ̃)sin(ψ̃) dφ̃=∫mconv(φ)ψ2g(ψ̃)sin(ψ̃) dψ̃,and the IVP becomes(43)m'conv(φ)=-f(φ)sin(φ)g(mconv(φ))sin(mconv(φ))mconv(φ1)=ψ2, φ1<φ<φ2.

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.

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 two-dimensional 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 piecewise-linear 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 O, 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)Ij=Pr(φj-1≤φ<φj)Δφ∫φ1φ2f(φ) dφ,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)Ij=Pr(φj-1≤φ<φj)Δφ∫φ1φ2f(φ)sin(φ)dφ.

More details regarding raytracing are available in [18], p. 34.

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)κ(θ):=-π-θ,θ<0,π-θ,θ>0.

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)ε(h,h*)=1N∑n=1N|hn-hn*|2,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 hrc*, respectively, where the “rc” subscript signifies “reconvolution”, that is hrc:=gdc*p.

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.

From the above list, it is clear that the step we propose in this paper, the deconvolution step, takes the least amount of time.

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 [α₁, α₂] when performing the (de-)convolution. We re-normalised p after truncation to ensure that ∫α1α2p(α) α=1. The problem is summarised in the box below, where(48)N(θ;μ,σ)=1σ2πexp-12θ-μσ2,represents the Gaussian, centred at μ with standard deviation σ. The value of φ₂ was chosen such that energy is conserved up to ∫φ1φ2fφdφ/∫ψ1ψ2g(ψ)d ψ≥10-3.

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.

The next step is to design the reflectors. In order to use the procedure outlined in Section 3, we need the limits ψ₁ and ψ₂. 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 piecewise-linear 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 [φ₁, φ₂]-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 Nr-1/2 behaviour of Monte Carlo raytracing ([18], p. 9), where N_r is the number of rays traced through the system.

There are a couple of regions where the raytraced distributions deviate from our predictions. Specifically, f* near φ1=-π/4 and φ2=29π/75, gdc* near κ(ψ2)≈1.26 and near both peaks of the Gaussians. The discrepancies in f* are due to the bins not aligning perfectly with the support of f, such that part of a collection bin may cross the φ₁- and φ₂-boundaries. The discrepancies of gdc* 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.

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 [ψ1,ψ2] and f is homogeneous on [φ1,φ2]. The scattering PDF p 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 h=g*p and attempt to recover g via deconvolution, then validate the reflectors we designed using gdc 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.

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² and 10⁴ 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⁻³ to 10⁻⁴ and 10⁻⁶. 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 worst-case scenario, as well as the best g_dc,, in the sense that it has the lowest RMS error.

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 (ψ), ψ ∈ [ψ₁, ψ₂], 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.

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⁴ iterations), respectively. In all these cases, we used a total of 10⁶ 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.

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 large-angle 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)L(θ;σ)=1πσσ2θ2+σ2,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 [α₁, α₂] and renormalised, such that ∫α1α2p(α)dα=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 g̃dc 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 h̃rc:=g̃dc*p with our prescribed h̃ in Figure 16. This may seem like a disappointing result, but let us compare h̃, h̃rc and h̃*p, 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 h 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 ε(h̃,h̃*p)=0.0624, whilst ε(h̃,h̃rc)=0.0356, i.e., our approach represents an improvement of approximately 40%. Visually, we see that the problematic regions for h̃rc are partly the peaks and partly near κ(γ1)≈1.02 and κ(γ2)≈2.22. 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.

We are now ready to design the reflectors. By fixing η = 0.001 and locating the two values of κ(ψ) where gdc(ψ)=η, we found κ(ψ1)=1.08 and κ(ψ2)=2.14, 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 three-dimensional version of the m_conv reflector. The raytraced distributions are shown in Figure 20a, where we see that the source sampling f̃* is correct, as is the resulting diffuse distribution h̃rc*. As for the intermediate specular target distribution, some deviations from the target g̃dc 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 105 rays. For the sake of completeness, we also raytraced the mdiv reflector in Figure 20b. Here, we only show g̃dc, g̃dc*, h̃rc and h̃rc*, 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 non-triviality 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.

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 in-plane 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 pre-processing 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 three-dimensional model of surface scattering using geometrical optics, with the ultimate goal of computing freeform reflectors with scattering surfaces. Lifting the requirement of in-plane 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.