An Improved Visible Normal Sampling
Routine for the Beckmann Distribution
Wenzel Jakob
August 2, 2014

Recently, Heitz and D’Eon [2014] proposed a method for importancesampling the distribution of visible normals in the context of microfacet BSDF models. One of their sampling routines internally relies
on a discontinuous mapping, which can cause problems in conjunction with Quasi Monte Carlo sampling and Markov Chain Monte
Carlo integration. In this report, we develop an alternative methodthat does not have this drawback.

1 Introduction
Microfacet BSDF models describe the interaction of light with random surfaces composed of microscopic dielectric or conducting facets that are oriented according to a
microfacet distribution. Validations against real-world measurements have shown thatmicrofacet models compare favorably against other families of parametric BRDF models, and for this reason they are a popular choice in the context of physically-basedrendering.
The standard approach [Walter et al. 2007] for importance sampling a microfacet
BSDF entails random sampling of a microgeometry normal proportional to a given
microfacet distribution, followed by the simulation of a reflection or refraction event.
However, due to masking effects, only a subset of the microgeometry receives lightfrom any particular direction. The standard technique does not take this into accountand generates many samples that contribute little or no energy, which can lead to highvariance in renderings.
Recently, Heitz and D’Eon [2014] proposed a method for importance sampling thedistribution of visible normals, which can lead to significant variance improvements,
particularly when dealing with very rough surfaces. However, one technical issue
still exists with their approach, specifically for surfaces modeled using the Beckmann

1

microfacet distribution: a key component of their algorithm maps a 2D uniform variate
to the x and y slopes of a surface facet. In their implementation, this mapping isdiscontinuous. While not a problem from a theoretical standpoint, this is undesirable
for two practical reasons:• It is common to use stratified samples or Quasi Monte Carlo integration to improve convergence in MC renderings. The effectiveness of these methods isgreatly diminished when scattering models use discontinuous mappings in their
the sampling routines.
• Some MCMC rendering techniques [Kelemen et al. 2002] expect that small perturbations to the input random variates will result in small changes to sampleddirections. This assumption is violated at discontinuities, which can lead to distracting visual artifacts due to variation in the speed of convergence in different
parts of the image.
The image below shows a plot of Heitz and D’Eon’s mapping from a 2D randomvariate to x facet slopes for the Beckmann distribution. The mapping was found using
a case-by-case analysis of the underlying problem, which gives rise to several “sheets”
with discontinuous transitions.

The goal of this report is the design of an alternative continuous mapping for thiscomponent.

1.1 Improved Sampling Routine
We only discuss sampling of the x slope (the y slope sampling does not need anymodification). The objective of this step is the inversion of the cumulative distributionfunction
G1 (ωi ) tan(θi )
1 12√
Cωi ( x ) :=exp − x + G1 (ωi )+ erf( x )2 22 π2on the interval [−∞, cot θi ], where G1 is Smith’s shadowing-masking term, ωi is theincident direction, and θi the associated latitude. For further details on this expression,refer to the supplementary material of the paper by Heitz and D’Eon [2014].
We first discard shared constant factors for simplicity:Ĉωi ( x ) :=tan(θi )√exp − x2 + erf( x ) + 1π

Let us switch to a new set of coordinates x 0 which relate to x as x 0 = erf( x ). In thesecoordinates, the above function is given bytan(θi )Ĉωi x0 = √exp −erf−1 ( x 0 )2 + x 0 + 1π
This reparameterization leads to a well-behaved function defined on a finite interval
x 0 ∈ [−1, erf(cot θi )]. Its shape varies based on the value of θi ∈ [0, π/2]. Two (renormalized and scaled) plots are shown below.

At perpendicular incidence, Ĉωi (x0 ) is a linear function (left), and for grazing θi it isslightly curved (right). Numerical root-finding techniques based on bracketing intervals can use its monotonicity property to reliably invert Ĉωi ( x 0 ) = ξ to solve for x 0given ξ ∈ [0, Cωi (erf(tan θi ))].
The derivative of this function is also simple to compute and given by
0Ĉω
x0 = 1 − tan(θi ) erf−1 ( x 0 )iDue the close-to-linear nature, Newton-type methods can be expected to converge ina low number of iterations. In Listing 1, we combine both approaches (bracketingintervals and Newton’s method) to safely invert Ĉωi given θi and a uniform variate ξ.
This listing also employs a heuristic that sets the starting point of the root-finding
iteration to a point that is close to the correct solution. This fit (obtained in Mathematica) ensures that only 1-2 iterations are required in practice.

3

void sample11(float thetaI, float U1, float U2, float &slope_x, float &slope_y) {const float SQRT_PI_INV = 1 / std::sqrt(M_PI);
float tanThetaI = std::tan(thetaI);
float cotThetaI = 1 / tanThetaI;/* Search interval (in the erf() domain) */float a = -1, c = erf(cotThetaI);/* Start with a good initial guess (approximation computed in Mathematica) */float fit = 1.0f + thetaI * (-0.876f + thetaI * (0.4265f - thetaI * 0.0594f));
float b = c - (1.0f+c) * std::pow(1.0f - U1, fit);/* Normalization factor for the CDF */float normalization = 1.0f / (1.0f + c +SQRT_PI_INV * tanThetaI * std::exp(-cotThetaI*cotThetaI));while (true) {/* Bisection criterion -- the oddly-looking boolean expression
are intentional to check for NaNs at little additional cost */if (!(b >= a && b <= c))b = 0.5f * (a + c);/* Evaluate the CDF and its derivative (i.e. the density function) */float invErf = erfinv(b);float value = normalization * (1.0f + b +SQRT_PI_INV * tanThetaI * std::exp(-invErf*invErf)) - U1;
if (std::abs(value) < 1e-5f)
break;/* Update bisection intervals */if (value > 0.0f)
c = b;else
a = b;/* Perform a Newton step */float derivative = (1 - invErf*tanThetaI) * normalization;b -= value / derivative;}/* Now convert back into a slope value */slope_x = erfinv(b);/* Y slope sampling works as before */slope_y = erfinv(2.0f * U2 - 1.0f);}

Listing 1: A C++ listing of the new sampling routine (single precision)

4

References
Heitz, E., and D’Eon, E. 2014. Importance sampling microfacet-based BSDFs using
the distribution of visible normals. Computer Graphics Forum (June).
Kelemen, C., Szirmay-Kalos, L., Antal, G., and Csonka, F. 2002. A simple androbust mutation strategy for the Metropolis light transport algorithm. Computer
Graphics Forum 21, 3, 531–540.
Walter, B., Marschner, S. R., Li, H., and Torrance, K. E. 2007. Microfacet models
for refraction through rough surfaces. In Proceedings of EGSR ’07, Eurographics
Association.

5