Rational BRDF

Romain Pacanowski, Oliver Salazar-Celis, Christophe Schlick, Xavier Granier,

Pierre Poulin, Cuyt Annie

To cite this version:

Romain Pacanowski, Oliver Salazar-Celis, Christophe Schlick, Xavier Granier, Pierre Poulin, et al..
Rational BRDF. IEEE Transactions on Visualization and Computer Graphics, 2012, 18 (11), pp.1824-
1835. �10.1109/TVCG.2012.73�. �hal-00678885�

HAL Id: hal-00678885

https://inria.hal.science/hal-00678885
Submitted on 14 Mar 2012

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 1

Rational BRDF
Romain Pacanowski, Oliver Salazar Celis, Christophe Schlick, Xavier Granier, Pierre Poulin and Annie Cuyt

Abstract—
Over the last two decades, much effort has been devoted to accurately measure Bidirectional Reflectance Distribution Functions (BRDFs) of
real-world materials and to use efficiently the resulting data for rendering. Because of their large size, it is difficult to use directly measured
BRDFs for real-time applications, and fitting the most sophisticated analytical BRDF models is still a complex task. In this paper, we introduce
Rational BRDF, a general-purpose and efficient representation for arbitrary BRDFs, based on Rational Functions (RFs). Using an adapted
parametrization we demonstrate how Rational BRDFs offer (1) a more compact and efficient representation using low-degree RFs, (2) an
accurate fitting of measured materials with guaranteed control of the residual error, and (3) an efficient importance sampling by applying the
same fitting process to determine the inverse of the Cumulative Distribution Function (CDF) generated from the BRDF for use in Monte-Carlo
rendering.

Index Terms—BRDF. Fitting. Importance Sampling. Monte-Carlo Rendering.

1 Motivation and Previous Work techniques. Of these, some focus on (hemi)spherical basis
functions, such as spherical harmonics [13], [14], Zernike
The BRDF is the keystone of the rendering equation [1], and polynomials [15], spherical wavelets [16], or spherical radial
thus of any lighting simulation process, since it describes the basis functions [17]. Others apply dimensionality-reduction
light reflection at a surface location. Consequently, building techniques such as homomorphic factorization [18], singular
a comprehensive BRDF representation suitable to explain, value decomposition [19], [20], or non-negative factoriza-
describe, and/or simulate all complex phenomena involved tion [21], [22]. Their main limitation resides in the fact that
in this process remains an active topic both in physics and the number of coefficients required to get an accurate result
computer graphics. The large number of scientific publications quickly grows with the directional frequency of the BRDF
dealing with BRDF can mainly be divided into two families, (quadratically, as shown by Mahajan et al. in [23]), and thus
whether they propose an analytical formulation, or rather a becomes intractable for highly specular BRDFs.
numerical process to approximate measured BDRF data.
Measured data may also be projected onto the parametric space
of analytical BRDF models. In fact, some models are specifi-
1.1 Closed-form vs. Numerical Representations
cally designed for numerical fitting [7], [24], [25]. Compared
to linear decomposition, fitting can easily handle materials
Papers in the first family present several models intended to with sharp specularity. However, this typically involves non-
represent some observed phenomena using either an empirical linear minimization techniques (e.g., the Levenberg-Marquardt
or a theoretical framework. These papers usually propose a algorithm), which are quite slow, and reaching a global
closed-form formulation for the BRDF, parametrized in a minimum of the function to be optimized is in general
manner that allows for the approximation of a limited set not guaranteed. Moreover, as analytical models have usually
of real-world reflectance behaviors. These models range from been designed for a specific class of materials, they cannot
efficient ad-hoc formulations such as Phong’s model [2], to accurately fit arbitrary BRDFs [26]. Using multiple lobes or
sophisticated ones that include the complex effects of wave combining different models may improve accuracy, but due
optics [3]–[5]. The most common group of models [6]–[11] to the many local minima, the optimization quickly becomes
exploit the micro-facet theory [12]. unpractical [21], [26].

Papers in the second family instead seek an efficient ap-

proach for representing the measured data using a set of 1.2 BRDF Importance Sampling
basis functions, by means of standard linear decomposition
Accurate fitting of measured data only solves half of the
• R. Pacanowski is with Univ. Bordeaux, LP2N, CNRS and IOGS.
problem. Plugging the fitted data into a rendering engine,
F-33400 Talence France. E-mail: romain.pacanowski@institutoptique.fr while preserving accurate and efficient computation is far from
straightforward. Over the years, Monte-Carlo-based techniques
• O. Salazar Celis and A. Cuyt are with the Department of Mathematics
and Computer Science, University of Antwerpen.
(see [27] for a comprehensive survey) have become the
standard approach for generating realistic images. As their
• C. Schlick and X. Granier are with the LaBRI - INRIA Bordeaux University. convergence speed is proportional to the inverse square root
E-mail: [schlick|granier]@labri.fr
of the number of samples and to the initial variance, it is
• P. Poulin is with the Department I.R.O., Université de Montréal. E-mail: critical to exploit efficient variance reduction techniques, the
poulin@iro.umontreal.ca most common one being Importance Sampling (IS) according
to the BRDF.
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 2

Original data size: BRDF = 99 MB, Tabulated CDF+PDF w30 MB Our approach: BRDF = 1.67 KB, Inverse CDF = 0.600 KB
Fig. 1. Monte-Carlo rendering with 2048 samples/pixel for a scene with three measured BRDFs from the MERL-MIT
database (blue-metallic on the dragon, beige-fabric on the floor, nickel on the sphere). Our approximation of the
BRDFs and the inverse CDFs, based on Rational Functions, provides efficient importance sampling with a negligible memory
footprint: with less than 1 KB of storage, our IS technique (right) offers equivalent quality (mean Lab difference is 0.77 and
max 7.03 on low-dynamic range images) compared to the reference solution (left) obtained from tabulated data of w 30 MB.
These tabulated CDF+PDF data have been generated by resampling the BRDF in (θv , θl , φl ) at 90 × 90 × 180. Furthermore,
the rendering time of our approach is 10% faster.

As detailed in Section 4.1, performing IS requires the in- are preferred in several numerical approximation problems in
verse of the BRDF’s Cumulative Distribution Function (CDF). scientific computing [31]. A Rational Function (RF) of a finite
There are basically two approaches to compute the inverse of dimensional vector x of real variables xi is:
the CDF. The first is to fit the measured data with an analytical n
X
BRDF model [7]–[10], [24], [28] that offers a readily invertible p j b j (xx)
CDF. In addition to the previously-mentioned weaknesses pn,m (xx) j=0
of non-linear fitting, all these approaches (except for [28]) rn,m (xx) = = m (1)
qn,m (xx) X
ignore the cosine factor that scales the BRDF according to qk bk (xx)
the incident light direction, reducing the efficiency of IS for k=0
grazing angles of light. The second approach consists of where the n+1 (resp. m+1) coefficients of the numerator (resp.
tabulating the CDF into a sorted data structure (e.g., binary denominator) are represented by the real numbers p j (resp. qk ),
search tree) and computing the inverse function on-the-fly and where b j (xx) and bk (xx) are multivariate basis functions.
in this structure [21], [29], [30]. A major benefit of this We use the multinomials in this paper, because they can be
approach is that the cosine scaling factor can be trivially evaluated efficiently. We order them by increasing total degree,
included, greatly improving the efficiency of IS. Unfortunately, for example in the bivariate case: b0 = 1, b1 = x1 , b2 = x2 , b3 =
the storage cost is several orders of magnitude higher than with x12 , b4 = x22 , b5 = x1 x2 , b6 = x13 , . . . Therefore, for a given degree
the first approach, and the iterative data retrieval process has we favor adding first smoother basis functions (e.g., x12 ) rather
a non-constant computation cost. than more oscillating ones (e.g., x1 x2 ). Furthermore, both
In this paper we introduce the following contributions: pn,m (xx)/qn,m (xx) and αpn,m (xx)/αqn,m (xx) take the same function
values for finite nonzero α, and the coefficients p j and qk need
• a general framework based on Rational Functions (RFs), only be determined up to a multiplicative constant that can
which efficiently represents BRDFs and CDFs without be used to normalize the representation of rn,m (xx). Therefore
having to separate diffuse and specular components. rn,m (xx) has no more than n + m + 1 free coefficients.
• an associated fitting technique that scales with the desired RFs are ideal for approximating data that exhibit steep changes
accuracy and memory footprint. The involved optimiza- which are characteristic for specular lobes. An illustration of
tion is that of a strictly convex function of which the approximation of lobe-like functions using RFs is given in
global minimum is guaranteed to be reached, provided Figure 2, where it can be observed that a low degree RF
that a feasible solution exists. can easily represent abrupt variations followed by regions
• a new Monte-Carlo estimator for importance sampling of almost constant values, whereas a polynomial with the
rendering, which does not require to store the PDF when same number of coefficients cannot. Such combinations of
combined with our representation. steep changes with flat regions are quite common in mea-
sured BRDF data and their corresponding CDF. However, in
computer graphics, RFs have seldom been employed (except
2 Rational Functions Framework for the ad hoc BRDF model proposed by Schlick [8]).
In approximation theory Rational Functions are recognized Algorithm 1 presents an overview of our fitting procedure
for their greater expressivity compared to polynomials. They based on the work of Salazar Celis et al. [32]. A preprocessing
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 3

Algorithm 1 Overview of the fitting procedure.

pdata = preProcessed(data)
Cmax = RF maximum number of coefficients
{Fi } = intervals for measured value x i of pdata
b0 , ..., bCmax −1 are ordered basis functions

1: fitData(pdata, Cmax , {Fi })

2: function fitData(data, Cmax , {Fi })

Fig. 2. Two low degree univariate approximations with 7 3: fitRatFunc(data, Cmax , {Fi })
4: while no solution do
unknowns: a polynomial of degree 6 (dashed line) and a 5: increase the width of all intervals {Fi }
rational function (full red line) consisting of a polynomial 6: solution = fitRatFunc(data, Cmax , {Fi })
of degree 1 and degree 5 in numerator and denominator 7: end while
respectively. The data are shown as dots. The RFs clearly 8: return solution
9: end function
follow the data much better and do not suffer from oscillations
2
like polynomials. Left: The rapidly decreasing function e−x , 10: function fitRatFunc(data, Cmax , {Fi })
sampled at 90 equidistant points. The maximum absolute 11: S = ∅ S best = ∅ ccn = ∞
error of the polynomial is 0.0689 while it is 0.00117 for the 12: for n = 0, ...,Cmax − 1 do
rational function. Right: The steeply decreasing lobe of the 13: for l = n, ..., 0 do
chrome-steel material from the MERL-MIT database. The 14: cn = ConditionNumber(A A(l, n − l))
15: S = Solve(P(l, n − l), data, {Fi })
maximum absolute error of the polynomial is 138 while it is 16: if S , ∅ and cn < ccn then
only 2.61 for the rational function. 17: ccn = cn
18: S best = S
19: end if
20: end for
21: end for
step is required before the fitting itself, for data reparametriza- 22: return S best
tion and noise removal. Contrary to classical non-linear meth- 23: end function
ods, our Rational Functions Framework (RFF) allows full
control of the residual error by setting intervals Fi around each
measured value and the involved optimization is guaranteed to problem P(n, m):
converge to a global minimum (cf. Section 2.1). Besides the
measurements, the memory budget (i.e., the maximum number arg min |cc|2
c ∈Rn+m+2
of coefficients), and the interval widths, the key part of the
subject to
algorithm requires solving a quadratic programming problem ( j)
P which is detailed in the following subsections. A(n,m
A n,m c − δ|A
j)
|2 > 0, j = 1, . . . , 2s + 2

with c = (p0 , . . . , pn , q0 , . . . , qm )t
( j)
and A n,m denoting the j-th row of matrix
 b0 (xx0 ) . . . bn (xx0 ) − f0 b0 (xx0 ) . . . − f0 bm (xx0 )
 
 . .. .. ..

2.1 Problem Statement for Fitting Data with RFs  .. . . .



 b0 (xx s ) . . . bn (xx s ) − f s b0 (xx s ) . . . − f s bm (xx s ) 

A n,m = 
 
−b0 (xx0 ) . . . −bn (xx0 ) f0 b0 (xx0 ) . . . f0 bm (xx0 ) 

Assume s + 1 measured values fi , each of them located at a  .. .. .. ..
 
 . . . .

vector x i . Let Fi = [ fi , fi ] represent a real-valued interval placed 
−b0 (xx s ) . . . −bn (xx s ) f s b0 (xx s ) . . .

around each fi . We would like to find a Rational Function f s bm (xx s )
rn,m (xx) that interpolates these intervals: and where | · |2 denotes the Euclidean norm. To avoid under-
or overflow of the computed coefficients, the real value δ > 0
is set to be the inverse of the conditional number of the matrix
∀i = 0, . . . , s fi 6 rn,m (xxi ) 6 fi , (2) A n,m . Furthermore, if P has a solution, it is unique and the
returned function rn,m is pole-free at each measured value,
i.e., ∀i = 0, . . . , s, qn,m (xxi ) > 0. Whether a solution exists for
a specific n and m, depends on the relation between the width
and this, for the smallest possible value for n+m, with n+m 6 s of the intervals Fi and the values n, m, and s. Algorithm 1 may
(usually n + m  s). find multiple solutions, but it always selects the more stable
solution by choosing the one associated with the lowest condi-
It can be shown [32] that a robust solution for rn,m (xxi ) can tion number of matrix A n,m . Although theoretically possible,
be obtained by solving the following quadratic programming multiple solutions were rare for BRDF/CDF fitting. Finally,
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 4

solving P can be done with any quadratic programming solver; Assuming that qn,m (xx) > 0, Inequality (3) becomes
in this paper, all RFs are obtained using the qpas routine of ∂pn,m pn,m (xx) ∂qn,m
the freely available QPC Matlab interface [33]. (xx) > (xx) .
∂x j qn,m (xx) ∂x j
In the next subsection, we introduce linear constraints that are By using Equation (2) the previous inequality is achieved for
easily added to P when fitting BRDFs or inverse CDFs. all sample positions x i if:
∂pn,m ∂qn,m ∂pn,m ∂qn,m
(xxi ) > (xxi ) fi and (xxi ) > fi (xxi ) . (4)
2.2 Additional Constraints
∂x j ∂x j ∂x j ∂x j
This last constraint consists of two additional linear inequal-
Value constraints: Since BRDFs are non-negative functions ities for the coefficients p0 , . . . , pn and q0 , . . . , qm per data
we build value constraints on the function as a whole into the point x i . The effect of adding these constraints can be seen
intervals right from the start: in Figure 3, where we show two fitted rational functions,
(left) without and (right) with the added discrete monotonicity
∀i = 0, . . . , s fi > 0 . constraints.

For 2D inverse CDFs (cf. Section 4) we make sure that:

1.5 1.5

∀i = 0, . . . , s 0 6 fi 6 fi 6 π/2 1 1

0.5 0.5
and similarly for 3D inverse CDFs:
0 0

∀i = 0, . . . , s 0 6 fi 6 fi 6 π .
1 1
0.8 1.5 0.8 1.5
0.6 0.6
0.4 1 0.4 1
0.2 0.5 0.2 0.5
µ 0 0 θν µ 0 0 θν
Even though this type of value constraints are only valid at
discrete positions, we found that it works satisfactorily. Fig. 3. Two low degree bivariate rational approximations
to CDFv−1 (µ) of the nickel MERL-MIT material. The ap-
Symmetry constraints: Helmholtz reciprocity (or symmetry)
proximation should be monotonically increasing in µ. Left:
is one of the main properties of physically based BRDFs
Rational interval interpolant obtained without Constraint (4),
and, by construction, inverse CDFs inherit it. Therefore it
clearly not monotonically increasing in µ. Right: Rational
is important to guarantee that the fitted RFs approximating
interval interpolant obtained with Constraint (4), now mono-
BRDFs or inverse CDFs preserve this property. Thanks to the
tonically increasing in µ.
parametrization that we are using for BRDFs (cf. Section 3.1)
we do not need additional constraints since the half-angle
parametrization already includes the symmetry property.
However, for inverse CDFs, we are using the classical light- 2.3 Solving P Efficiently
view parametrization (cf. Section 4.1), and therefore we
The speed for the resolution of P greatly depends on the
enforce symmetry by grouping basis functions that should
number s + 1 of intervals to be interpolated. Furthermore, the
receive the same coefficients. For example, assume we would
faster P is solved, the sooner the algorithm is going to find the
like to fit data with a trivariate RF:
optimal solution and return. To quickly find the couple (n, m)
rn,m (x1 , x2 , x3 ) = c0 x13 x24 x3 + c1 x23 x14 x3 (c0 , c1 ) ∈ R2 we propose an adaptive procedure that operates on a selected
small subset of the whole data. The key observation here is
and enforce its symmetry for the variables x1 and x2 , i.e., that if the current tested rn,m (xx) RF is not a solution for the
rn,m (x1 , x2 , x3 ) = rn,m (x2 , x1 , x3 ). Instead of adding, into P, selected subset, it cannot be for the whole dataset.
equality constraints between c0 and c1 , we group the two basis
functions into one: Of the given s + 1 data intervals, only a small number s0
 of intervals are uniformly selected and rn(0) 0 ,m0
(xx) is computed


rn,m (x1 , x2 , x3 ) = c0 x13 x24 + x23 x14 x3 . such that it satisfies Equation (2) for these s0 data; these s0
intervals are called the training data. Then it is checked how
This reduces the size of matrix A n,m and hence the size of the many of the original s + 1 interval interpolation conditions
quadratic programming problem. are automatically satisfied by rn(0) 0 ,m0
(xx) in addition to the s0
imposed conditions. Usually this is quite a lot more. These
Monotonicity constraints: Since inverse CDFs are positive
s + 1 − s0 intervals are called the verification data. Among
monotonic functions it is crucial to enforce that the derivative
the violated interval interpolation conditions, we select s1 − s0
of the fitted function is non-negative. In other words, we want
(s1 − s0 = 1 in our implementation) additional data points to
to ensure the monotonicity of the rational function rn,m (xx) =
compute rn(1)
1 ,m1
(xx) that satisfies Equation (2) for these s1 data.
pn,m (xx)/qn,m (xx) in one of the variables x j :
In other words, we update the set of training data. These s1 − s0
∂rn,m additional training data are selected where rn(0) 0 ,m0
(xxi ) deviates
(xx) > 0. (3)
∂x j most from the given intervals Fi . With rn(1) 1 ,m1
x
(x ) we then check
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 5

the new s + 1 − s1 verification data again. And so on, until the

obtained rational model satisfies all verification data.

3 Rational Functions for BRDF

3.1 Parametrization
Fig. 4. Lens-flare-like acquisition artefacts of the grease-
A BRDF ρ is a four-dimensional non-negative function defined covered-steel MERL-MIT material illuminated by 9 direc-
on the Euclidean product of two hemispheres and depending tional sources. Left: Original 3D data. Middle: 3D data
both on the lighting l and viewing v directions. The most without noisy grazing angles (i.e., θh or θd > 80◦ ). Right: 2D
common way to express l and v is to use the standard zenith projection on (θh , θd ) mostly removes all artefacts.
and azimuthal angles (θl , φl ) and (θv , φv ): ρ(θl , φl , θv , φv ), often
denoted by ρ(ll,vv).
Although most acquisition devices (e.g., gonioreflectometers) As illustrated in the next sections, a single low-complexity
cover the acquisition space by uniformly sampling these RF is able to approximate isotropic BRDF data well, with-
coordinates, the computer graphics community prefers to use out requiring multiple lobe fittings like the current standard
the BRDF parametrization introduced by Rusinkiewicz [34]: techniques, even when several phenomena (diffuse reflection,
ρ(ll,vv) = ρ(hh,dd ) = ρ(θh , φh , θd , φd ) where h is the half-vector forward or backward glossy or specular reflections, Fresnel
used in the micro-facet theory, and (θd , φd ) are the spherical effects, etc.) are simultaneously observed on a given material.
coordinates of vector l expressed in a rotated frame where h Although isotropic BRDF measurements have been intensively
defines the north pole direction. This parametrization offers studied, much less work has been done on fitting measured
several interesting benefits: (1) in the case of isotropic ma- anisotropic materials, basically because much less anisotropic
terials (which are invariant when rotating the surface around data is publicly available. In Section 5 we present our adap-
its normal vector), one can set φh = 0 without any loss of tation of Rational BRDF to an anisotropic data model, and
generality, (2) the reciprocity of light propagation ensures discuss the need for high quality measured 4D datasets.
that φd can be limited to the range [0, π], again without any
loss of generality, (3) as noted by Romeiro et al. [35], for a
very common class of materials called bilateral symmetric, the 3.2 Rational Approximation of BRDF
domain of φd can even be reduced to the range [0, π/2], and
finally (4) as several phenomena involved in light reflection The most comprehensive and accurate database of measured
are mostly decorrelated along the different parameter axes, it isotropic BRDFs, known as the MERL-MIT BRDF database,
is relatively safe to factor the 4D function into a product (or combines over 1 billion individual BRDF measurements gen-
a sum of products) of two 2D functions, as done in work on erated by Matusik et al. [20]. The measured BRDF data are
factorization of BRDFs and SVBRDFs [21], [36]. available with a 90 × 90 × 180 angular sampling in (θh , θd , φd ),
which represents a storage amount of 33 MB per material
As the BRDF is intrinsically 3D/4D for isotropic/anisotropic
(the database was last updated in 2006). We have tested
materials, further dimension reduction cannot be achieved
our RF approximation technique on many different materials
without generating some undesired folding of its domain.
from the database, but for the clarity of presentation, we
This can be seen as a standard aliasing process, where an
focus here on four representatives of common BRDF fam-
infinite number of configurations of the 4D space are projected
ilies: beige-fabric (almost perfect Lambertian reflection),
onto a single configuration in the lower dimensional space.
blue-metallic-paint (glossy reflection with strong chro-
However some recent work has shown that for most isotropic
matic behavior), nickel (specular reflection), chrome (almost
materials, the BRDF can be projected onto a well-chosen 2D
perfect mirror reflection), which are presented in Figure 7 (top
parametric space without severe visual degradation [35], [37].
row).
Romeiro et al. [35] have observed that almost all materials of
the MERL-MIT database [20] are visually well approximated For each of these materials, we fit, as described in Section 2,
by a projection onto (θh , θd ) (see Figure 7, top and middle a bivariate rational function rn,m (θh , θd ) that satisfies Equa-
rows). Incidentally, we have observed that the filtering operator tion (2). Of course, the behavior of the approximation is de-
involved when projecting measured data on (θh , θd ) sometimes
h i
termined by the choice of the intervals Fi = fi , fi containing
also results in more visually pleasing rendering results (see ρ(θhi , θdi ), (i = 0, . . . , 90 × 90 − 1), for each of the BRDF values.
Figure 4). Ideally, the widths of these intervals are calibrated such that
In this paper, we propose to go a step further and show that they respect the accuracy of each individual measurement.
we can accurately approximate the whole 2D projection of Since measurement error bounds are not explicitly available
the BRDF ρ(θh , θd ) by a single RF rm,n (θh , θd ), defined as in here, the interval widths are taken such that the renderings
Equation (1), that we call a Rational BRDF: are visually satisfactory, while leaving a reasonable number
m + n of coefficients (cf. Figure 6) in the approximation. As
ρ(ll,vv) ≈ ρ(θh , θd ) ≈ rn,m (θh , θd ). (5) a general guideline, we control the interval widths with the
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 6

blue-metallic-paint terials. The environment map is approximated by using 1024

1.2 0.35 directional light sources selected according to their power in
1
0.3
the surrounding environment map. As can be observed, the
0.8
0.25
visual error (measured in Lab color space) is approximately
0.2

0.6 of the same magnitude as with the 2D projected data, but

0.15

0.4
0.1
our Rational BRDFs require between 165 and 1000 times
0.2
0.05
less storage memory (in our current implementation, each
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4
coefficient p j and qk of the RF is stored as an 8-byte double-
θd (θh = 0) θh (θd = 0) precision float). When compared to the original 3D data,
nickel the compression rate ranges from 28700:1 to 174000:1. Note
400
that using 4-byte standard floats is likely to be sufficient in
5

350
many cases. Finally, for grazing angles configuration and for
300
4
rendering purposes only, we extend our RF by clamping light
250
3 and view directions that exceed 75 degrees. As shown in
200

150 2
Figure 1 this does not introduce visual artefacts. For all our
100
figures, we used the same clamping for the data and the RFs
1
50 in order to illustrate the accuracy of the fitting.
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4

θd (θh = 0) θh (θd = 0) The final approximation errors are summarized in the table
chrome of Figure 8: we have computed for (column 1) each material
350
(column 2) the maximum relative error of the 2D data using
2500

300
the (θh , θd ) parametrization, our 2D RF fitting against (col-
2000
250 umn 3) the 2D data and against (column 4) the original 3D
1500 200 data, and (column 5) using a non-linear fitting procedure with
1000
150 the modified Ashikhmin-Shirley BRDF model as described
100
by Ngan et al. [26]. As expected, for every material and
500
50
regardless of the representations, the fitting errors grow with
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4
the specularity of the material. As can be observed, the overall
θd (θh = 0) θh (θd = 0)
errors of the rational approximations and the Ngan fitting are
Fig. 5. Specular materials’ BRDF data (dots) without grazing quite similar. However the interval interpolation process has
angles (i.e., θh or θd > 75◦ ) and corresponding RF approxima- the advantage to offer full control both on the location and
tions (full line) per color channel (RGB). the magnitude of the errors. We chose to compare with the
Ashikhmin-Shirley model because, as stated by Ngan et al.
following formula: [26]: “for single specular lobe, Cook-Torrance, Ashikhmin-
Shirley and the He models perform well for most of the
Fi = [ fi , fi ] = [ρ(θhi , θdi )(1 − i ), ρ(θhi , θdi ) (1 + i )] 100 isotropic materials.” Furthermore, the effect of adding
where i controls the desired relative error. Figure 8 presents one lobe (hence fitting with a two-lobe BRDF model), as
the maximum relative errors obtained for BRDFs. However, shown by Ngan et al. [26] (in their Section 4.2), reduces
the interval widths (and therefore the i ) are most of the time indeed the fitting error (by approximately 25%) for 26 of the
not set uniformly. This is basically done by choosing smaller 31 materials. However, according to these authors, the fitting
interval widths near the hemisphere pole and more relaxed process quickly becomes unstable with more than three lobes
widths near grazing angles. and the benefit of doing so is only marginal. In conclusion, the
errors presented in the table of Figure 8 for the Ashikhmin-
For all materials we noted increased acquisition noise when Shirley model has been improved at most by 25% when adding
approaching grazing angles (θh > 60◦ and/or θd > 60◦ ). In a second lobe. Therefore, a model from a two-lobe fitting
the case of highly specular materials, we also noted that the would still not be competitive (regarding the error and the
acquired BRDF values around the center of the lobe (θh = 0) numerical stability of the technique) against our RF fitting
are particularly large and subject to noise as illustrated in technique.
Figure 5 on chrome and nickel BRDFs. Small interval
Finally, the rightmost column (column 6) of the table of
widths are thus chosen on data with low noise, leading to
Figure 8 shows the maximum relative error against the 3D
a very accurate fit in such a configuration (as also shown in
data when using least-squares bivariate polynomials with
Figure 5). However, as confirmed by Figures 1 and 7, even in
the same number of coefficients as the RFs. Except for
the case of materials including very noisy data and therefore
the low frequency beige-fabric material, the errors are
larger interval widths, the rendering obtained with the fitted
several orders of magnitude higher than the ones obtained
function has a very low visual impact when compared to the
with RFs. Furthermore, these errors show that approximating
rendering based on the original data.
high frequency signals with RFs is more efficient than using
Figure 7 (bottom row) presents the rendering obtained with polynomials. As shown by Mahajan et al. [23], the required
the Rational BRDF approximation of our four selected ma- number of coefficients for polynomial functions (e.g., spherical
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 7

harmonics) grows quadratically with the directional frequency PDFv (θl ) CDFv−1 (µ)
90
of the BRDF. Note that it is not straightforward to fit directly 100 80
spherical harmonics (SHs) using the (θh , θd ) because these 70
80
angles do not span a spherical domain. Other approaches 60
50

pdf
such as Sillion et al. [38] or Westin et al. [14] have used 60

θl
40
classical spherical parametrization but require between 80 and 40 30

133 coefficients for diffuse or glossy materials, and become 20

20
10
impractical for high frequency BRDF. 0 0
0 20 40 60 80 0 0.2 0.4 0.6 0.8 1
θl µ
Material (n + 1, m + 1) Maximum degree
Fig. 9. Comparisons of the PDFv and the inverse CDFv−1 for
R G B R G B
nickel material. The measured inverse CDF curves traced
beige-fabric (18,12) (21,12) (18,12) (5,4) (5,4) (5,4)
in red for θv = 0◦ , in green for θv = 45◦ , and in blue for θv = 85◦ ,
blue-metallic (20,28) (22,24) (16,16) (5,6) (6,6) (5,5)
are always simpler than the ones for the PDF.
nickel (27,16) (40,9) (36,12) (6,5) (8,3) (7,4)
chrome (84,17) (47,14) (50,26) (12,5) (9,4) (10,6)
Fig. 6. The data for the fitted RF representation for each Starting from a measured BRDF scaled by n ·ll, in other words
material 2D BRDF are provided for the denominator and the ρ(ll,vv) cos θl , we first compute a tabulated version of the inverse
numerator, each per (R,G,B) color channel: (center column) CDFs:
number of coefficients, (right column) maximum degree. • a 2D table for CDF−1 (θv , µ) where θv ∈ [0, π/2] indexes
the view direction and µ ∈ [0, 1] is a random sample,
• a 3D table for CDF−1 (θv , τ|θl ) where θl ∈ [0, π/2] and τ ∈
4 Rational Functions for CDF [0, 1] is a second random sample.

4.1 Parametrization Since we are using isotropic materials there is no azimuthal

dependency on the view direction v = (θv , φv ), i.e., φv = 0. Then
As mentioned in Section 1.2, efficient Importance Sampling we approximate these tabulated CDF data with a bivariate and
(IS) strategies significantly increase the rate of convergence trivariate RF respectively:
of Monte-Carlo rendering techniques. In this section we
 θl = CDF (θv , µ) ≈ rnθ ,mθ (θv , µ)
 −1


present quasi-optimal IS of arbitrary BRDFs by means of
 φl = CDF−1 (θv , τ|θl ) ≈ rnφ ,mφ (θv , θl , τ) .


RF approximation. With brevity in mind, we only detail our
procedure for the standard (ll,vv) parametrization. Note that this Note that since the CDF is the integral of its corresponding
parametrization has been proven to be well-suited for glossy PDF, the shapes of the inverse CDFs are always much simpler
surfaces [21]. than the PDFs, as shown in Figure 9, and RFs can easily
The principle of IS rendering is to define a stochastic estimator reproduce the steep variations that are present in most datasets.
for the reflected radiance L(vv) in direction v by averaging The gain of this approach is twofold. First, due to the compact
the contribution from K random light directions l k selected representation of RFs, generating each random vector l k is very
according to a conditional Probability Density Function (PDF) efficient as it only involves the evaluation of low-degree mul-
PDF(llk |vv) = PDFv (llk ): tivariate RFs. This is actually done in constant time, contrary
K to approaches based on the numerical inversion of tabulated
1 X n ·llk
L(vv) ≈ ρ(llk ,vv) L(llk ). (6) CDFs. Second, it allows for a more effective IS strategy
K k=1 PDFv (llk ) by using our new Monte-Carlo estimator (the derivation is
available in the dedicated companion document, Section 2),
When there is no prior knowledge about the incident lighting, which is directly based on the inverse CDF instead of the
the optimal choice for the estimator is to use a PDFv (llk ) that PDF:
is proportional to the BRDF scaled by n · l k . Each random K
1 X
light direction l = (θl , φl ) required by the estimator can then L(vv) ≈ αv (µk , τk ) ρ(vv,llk ) (nn ·llk ) L(llk )
be obtained by generating a pair of uniform random numbers K k=1
(µ, τ) and inverting a pair of Cumulative Distribution Functions ∂CDFv−1 ∂CDFv−1 (7)
(CDFs), obtained from the integration of the selected PDF: with αv (µ, τ) = (µ) (τ | θl ) sin θl
∂µ ∂τ
θl = CDFv−1 (µ) and φl = CDFv−1 (τ|θl ). and θl = CDFv−1 (µ) .
Since both inverse CDFs are RFs, their partial derivatives can
We propose an innovative technique intended to combine
easily be evaluated on the fly.
the strengths of previous approaches. The idea is to directly
define a closed-form expression for the inverse CDF, without Unlike previous approaches (e.g., [21]), which use different
preliminary analytical formulation of either the CDF or the parametrizations depending on the type of material, we have
PDF. decided to use the (θl , φl ) parametrization for all materials. This
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 8

beige-fabric blue-metallic-paint nickel chrome

Data size: 33 MB Data size: 33 MB Data size: 33 MB Data size: 33 MB

Data size: 190 KB [183:1] Data size: 190 KB [183:1] Data size: 190 KB [183:1] Data size: 190 KB [183:1]
Max Lab error: 4.65 Max Lab error: 5.13 Max Lab error: 17.1 Max Lab error: 19.7
Mean Lab error: 1.00 Mean Lab error: 1.33 Mean Lab error: 2.85 Mean Lab error: 1.19

Data size: 0.85 KB [40500:1] Data size: 1.15 KB [34170:1] Data size: 0.19 KB [91125:1] Data size: 0.33 KB [52071:1]
Max Lab error: 9.07 Max Lab error: 6.51 Max Lab error: 18.0 Max Lab error: 46.5
Mean Lab error: 1.23 Mean Lab error: 1.38 Mean Lab error: 1.52 Mean Lab error: 1.90
Fig. 7. From top to bottom: original measured data in (θh , θd , φd ), projected data in (θh , θd ), rational BRDF fit in (θh , θd ). All
scenes are rendered using 1024 directional light sources sampled from the surrounding environment map. Each Rational
BRDF generates a similar magnitude of Lab visual errors as the ones produced by the 2D projected data, but requires
between 165 and 1000 times less storage.

provides a consistent choice and we leave for future work the a 91 × 91 grid. Similarly, the values of CDFv−1 (τ|θl ) ∈ [0, π] are
question of the most suitable parametrization for importance computed on a 91 × 91 × 128 grid, which size is then reduced
sampling independent of the material type. Nevertheless, as by a factor of two using the fact that the inverse CDF is
illustrated in Section 4.2, the whole IS strategy proposed symmetric. Second, to accelerate the resolution of P, we use
here results in a very efficient computation, whatever the 300 measured values as training data (cf. Section 2.2).
complexity of the underlying BRDF, with extremely compact
storage compared to prior work.

As explained in Section 2.2 we guarantee that the fitted RF for

4.2 Rational Approximation of CDFs the 3D inverse CDF is symmetric with respect to θv and θl ,
i.e., rnφ ,mφ (θv , θl , τ) = rnφ ,mφ (θl , θv , τ). Moreover, we guarantee
We proceed as follows. First, we compute two tabulated that both RFs rnθ ,mθ (θv , µ) and rnφ ,mφ (θv , θl , τ) are monotonically
versions of the inverse CDFs for each of the four selected increasing with respect to µ and τ, respectively. Besides mono-
materials. The values of CDFv−1 (µ) ∈ [0, π/2] are computed on tonicity and symmetry properties, we also enforce boundary
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 9

2D proj. Romeiro Rational Functions Non-linear A&S Polynomial Functions

33 000 KB 0.19 KB – 1.15 KB 0.0625 KB 0.19 KB – 1.15 KB
Error on 3D data Error on 2D data Error on 3D data Error on 3D data Error on 3D data
beige-fabric (0.71, 0.73, 0.69) (0.71, 0.60, 0.55) (0.71, 0.71, 0.69) (1.66, 1.75, 1.68) (0.71, 0.72, 0.71)
blue-metallic (0.80, 0.80, 1.0) (0.25, 0.42, 0.29) (0.74, 0.71, 1.01) (28.49, 34.80, 17.50) (43.47, 77.77, 69.02)
nickel (1.12, 2.44, 2.95) (0.77, 0.78, 0.86) (0.77, 2.45, 1.95) (4.99, 4.11, 4.64) (47 697, 49 399, 152 752)
chrome (5.67, 4.01, 3.99 ) (1.72, 1.97, 3.36) (3.29, 2.38, 3.38) (65.04, 133.63, 148.42) (26 941, 1 323 656, 1 125 357)
Fig. 8. Maximum relative error of the approximation, per (R,G,B) color channel. In the usual order: 2D projected data
compared to 3D original data; rational BRDF compared to the 2D projected data, and compared to the 3D original data; non-
linear fitting of the Ashikhmin-Shirley analytic model compared to the original data, and least-squares bivariate polynomials
compared to the original 3D data. All data at grazing angles (> 80◦ ) have been discarded.

conditions of the CDFs: where rn,m is the corresponding rational approximation. Note
from these figures that the number of coefficients needed is
CDF (θv , µ = 0) = 0
−1
CDF−1 (θv = 0, θl , τ) = πτ
π always very moderate. This indicates that the proposed forms
CDF−1 (θv , µ = 1) = CDF−1 (θv , θl = 0, τ) = πτ in Equations (8) and (9) are adequate and capture most of
2
the information. The table of Figure 14 presents the fitting
by imposing the following rational forms: time for each of the 2D CDFv−1 (µ). These timings are directly
π pn ,m (θv , µ) proportional to the number of coefficients tested.
rnθ ,mθ (θv , µ) = µ + µ(1 − µ) θ θ (8)
2 qnθ ,mθ (θv , µ)
Figure 13 further compares, for each material, the rendering
pnφ ,mφ (θv , θl , τ) quality obtained (first row) with tabulated CDFs and PDFs
rnφ ,mφ (θv , θl , τ) = π τ + τ(1 − τ) θv θl . (9)
qnφ ,mφ (θv , θl , τ) against (second row) approximated CDFs. We also compute
the variance of the sphere pixels in each image to compare
For each of the four selected materials, we apply the fitting the performance of the IS strategy based on our RF CDF with
algorithm described in Section 2, complemented with the that of the tabulated data. To provide a fair comparison of the
corresponding modifications and new constraintsh mentioned i sampling efficiency, our renderer uses the BRDF combined
in Section 2.2. The choice of the intervals Fi = fi , fi , now either with the tabulated CDF or its RF approximation. For the
containing either θli or φil , for each of the inverse CDF values, blue-metallic, beige-fabric, and nickel materials the
completely determines the behavior of the approximations. resulting variance is as low as the reference solution obtained
Since these inverse CDFs guide the IS process, the interval with tabulated data. For the chrome material the variance is
widths are taken such that the resulting variance is acceptable, 10× larger than the reference solution, but the Lab difference
while leaving the number m + n of coefficients in the approxi- between the images shows that our RF approximation still pro-
mation manageable. As illustrated below, this basically implies vides a visually satisfactory result. Chrome remains the most
a fairly accurate CDFv−1 (µ) approximation, and capturing the challenging material. However, to our knowledge, previous
main trend of CDFv−1 (τ|θl ). For the 2D CDFv−1 (µ) we define approaches could not even represent the chrome data.
the intervals as:
Note that our technique guarantees the existence of a rep-
 fi = θl − (1 + θl )
i i

 resentable inverse CDF without introducing a bias in the
i = 0, . . . , 91 × 91 − 1

 fi = θi + (1 + θi )

 sampling process because the RFs preserve monotonicity.
l l
However, the speed of convergence of the IS is directly
with  = 0.015 for beige-fabric,  = 0.035 for dependent on the overall quality of the fitted function.
blue-metallic,  = 0.05 for nickel, and  = 0.1 for
the highly specular chrome. In a similar way, we define the Regarding the size of our functions, Figure 10 shows that
intervals for 3D CDFv−1 (τ|θl ) according to φl with  = 0.015 for with RF approximations, we achieve high compression rates.
beige-fabric and blue-metallic,  = 0.72 for nickel, For example, the nickel material in our approach is 7600×
and  = 0.8 for chrome. For the 3D CDFv−1 (τ|θl ), we relax (resp. 3100×) more compact than the one proposed by [21]
the intervals near θv = θl and near grazing angles θv , θl > 80◦ . (resp. [30]). The numbers from the first row of the table
In order to keep the number of coefficients manageable in have been directly reported from the original paper whereas
rnφ ,mφ (θv , θl , τ), we restrict the partial degree in τ to at most for the other rows, we have reported the lowest available
2, which is sufficient for capturing the steep increase in τ. sizes, which are coming from Montes’ Ph.D. thesis [39].
Notice that previous approaches do not provide any results for
Figure 11 shows the 2D approximations of CDFv−1 (µ) whereas
specular chrome whereas our technique may still be used to
Figure 12 shows the 3D approximations of CDFv−1 (τ|θl ). The
approximate an importance function. This high compression
relative errors shown are
ratio directly results from the properties of RFs, which can
|CDF−1 − rn,m | efficiently represent steep changes with low-degree approxi-
1 + |CDF−1 | mations, as well as from the fact that we do not store any
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 10

beige-fabric blue-metallic nickel chrome

Factored na ' 200 ' 200 na
Cascade PDF 374.16 1685.64 2057.51 na
A. Disc 272.62 715.91 846.15 na
RFs 0.117 0.211 0.230 0.172
Fig. 10. Comparison of the CDF and PDF data size in KB for
different techniques: Factored representation [21], Cascade
PDF [29], A. Disc [30], and Rational Functions. Our Rational Fig. 15. A blue-metallic-paint sphere illuminated by
Function approach consumes always less memory than the a single directional light source, with different anisotropic
previous approaches. transformations using a univariate function rma0 ,n0 (φh ). The
original isotropic material is shown on the left. The middle
beige-fabric blue-metallic nickel chrome
image corresponds to an elliptical scaling and presents the
2.82 5.74 19.19 6.44 typical butterfly pattern of brushed metals. The image on the
Fig. 14. Fitting timings in seconds computed with our Matlab right corresponds to a star-shaped scaling function.
program on a Intel Xeon L5420@2.5GHz for the 2D inverse
CDF. These timings include all iterations of the fitting algo-
rithm detailed in Section 2. results obtained with different anisotropic transformations of
the blue-metallic-paint material.
In theory, there is nothing preventing our RF method from
PDFs. As explained earlier, the PDF is computed on the fly approximating a measured full 4D BRDF. The results that
from the inverse CDF during the rendering step and thus does we obtained so far suggest that the 3D subset (θh , θd , φh )
not require any supplemental storage. of Rusinkiewicz’s parametrization is sufficient for some
Finally, the efficiency of the IS process can be observed in anisotropic BRDFs, similar to the use of (θh , θd ) for isotropic
Figure 1 which presents a global illumination rendering of BRDFs [35]. Reducing the dimension to 3D leads to less
a more complex scene. Compared to the reference solution, constraints for the interval interpolation technique [32], and
resulting in a small Lab difference, our solution achieves the results in a large speed-up of the fitting process. Unfortunately,
same overall rendering quality but with a drastically lower far less data is publicly available. We experimented with
memory footprint (90 MB + 30 MB vs. 1.67 KB + 0.6 KB) another MIT BRDF database [26] containing four measured
when using the 3D measured data, and (570 KB + 30 MB vs. anisotropic materials. However as noted in previous work [10],
1.67 KB + 0.6 KB) when using the 2D projected data. the embedded data contain significant acquisition noise, and
present large areas in its parametric domain where the mea-
sured BRDF values are quite sparse. These two issues lead to
5 Anisotropic BRDF Model using RFs artifacts in the fitted data. Smoothing and good quality data
completion [11] should improve our results, but this remains
Ward [7] observed that for many anisotropic materials, the one of our future research directions.
variation of the reflected intensity when rotating the surface
around its normal vector often consists of a simple scaling
factor applied to an average isotropic lobe. Based on that 6 Conclusion and Future Work
observation, we propose the following model for anisotropic In this paper, we have introduced Rational BRDFs, a generic
BRDFs as a product of two rational functions: and compact representation of arbitrary BRDFs, based on
ρ(l, v) ≈ rma0 ,n0 (φh ) rm,n
i
(θh , θd ) Rational Functions. By employing a subset of the BRDF
parametrization introduced by Rusinkiewicz [34], isotropic
where rm,n i (θ , θ ) represents a standard isotropic Rational
h d BRDFs can be represented as bivariate RFs. Very compact
BRDF and rma0 ,n0 (φh ) is a scaling factor defining the anisotropic approximations of measured BRDFs are obtained with a mem-
variation. This scaling factor agrees with classical brushed- ory footprint that is usually less than one kilobyte for arbitrary
metal behavior [7] and with several other anisotropic mod- isotropic materials. Moreover, the same approximation process
els [8], [10], [28]. applies to the inverse CDF generated from the corresponding
BRDF multiplied by the cosine factor. Combined with our new
To validate this model, we have chosen to generate “hybrid”
Monte-Carlo estimator which exploits the RF formulation, our
anisotropic BRDFs by post-processing measured isotropic
approach offers a quasi-optimal IS scheme, with compact stor-
materials. We may view the scaling function as an approxima-
age compared to prior work. We have shown that the Rational
tion of the directional variation generated by weaving meso-
BRDF is suitable to reproduce some anisotropic effects on a
structures of textiles, or regular geometric structures of crys-
3D subspace of Rusinkiewicz’s parametrization [34] by the
tals. Actually, the whole process can be seen as a combination
introduction of a univariate rational function that scales an
of the hybridization technique, based on a linear combination
isotropic Rational BRDF according to its orientation.
of materials, proposed by Matusik et al. [20], and the meso-
structure simulation technique used by Westin et al. [14] to All of our results are the foundations for future research. First,
predict reflectance from complex materials. Figure 15 presents we have started investigating acquisition techniques to obtain
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 11

beige-fabric blue-metallic-paint nickel chrome

Mean relative error: 0.0035 Mean relative error: 0.012 Mean relative error: 0.0105 Mean relative error: 0.0343
Max relative error: 0.0119 Max relative error: 0.0305 Max relative error: 0.0416 Max relative error: 0.0978
Nb coefficients: (6,8) Nb coefficients: (11,8) Nb coefficients: (18,24) Nb coefficients: (11,16)
Fig. 11. Rational Function approximation of the bivariate CDFv−1 (µ) for our four selected MERL-MIT materials. For each
material, the colormap represents the fitting relative error and ranges from the minimum relative error (black) to the maximum
relative error (yellow). The number of coefficients is expressed as a pair (numerator,denominator).
beige-fabric blue-metallic-paint nickel chrome

Mean relative error: 0.020 Mean relative error: 0.030 Mean relative error: 0.070 Mean relative error: 0.099
Max relative error: 0.137 Max relative error: 0.138 Max relative error: 0.713 Max relative error: 0.774
Nb coefficients: (9,7) Nb coefficients: (19,21) Nb coefficients: (7,10) Nb coefficients: (7,10)
Fig. 12. Slices (θv = 45◦ ) of the Rational Function approximation of the trivariate CDFv−1 (τ|θl ) for our four selected MERL-MIT
materials. For each material, the colormap represents the fitting relative error and ranges from the minimum relative error
(black) to the maximum relative error (yellow). The mean and max relative errors are given for the whole 3D data. The
number of coefficients is expressed as a pair (numerator,denominator).

dense BRDF measurements of materials, focusing on families Acknowledgments

of BRDFs that are currently lacking in public databases such
as materials exhibiting retro-reflective behavior, and various References
types of anisotropic reflections. Based on this dataset, we may
extend the previous studies on parametrization [35], [37] to [1] J. T. Kajiya, “The rendering equation,” in Proc. ACM SIGGRAPH ’86,
more generic materials. In investigating extensions to the work 1986, pp. 143–150.
of Lawrence et al. [21], we also hope to find a set of best-
[2] B. Phong, “Illumination for computer generated pictures,” Comm. ACM,
suited for IS. vol. 18, no. 6, pp. 311–317, 1975.

[3] X. He, K. Torrance, F. Sillion, and D. Greenberg, “A comprehensive

physical model for light reflection,” in Proc. ACM SIGGRAPH ’91, 1991,
pp. 175–186.

Second, as shown by the anisotropic transformation of [4] J. Stam, “Diffraction shader,” in Proc. ACM SIGGRAPH ’99, 1999, pp.
75–84.
blue-metallic-paint presented above, generating hybrid
materials by combining analytical and measured BRDFs ap- [5] I. Icart and D. Arquès, “A physically-based BRDF model for multilayer
systems with uncorrelated rough boundaries,” in Proc. EGWR ’00, 2000,
pears to be a promising direction of research. We are currently pp. 353–364.
formalizing this process to define a flexible framework for
intuitive BRDF editing. [6] R. Cook and K. Torrance, “A reflectance model for computer graphics,”
ACM Trans. Graph., vol. 1, no. 1, pp. 7–24, 1982.

[7] G. Ward, “Measuring and modeling anisotropic reflection,” in Proc.

ACM SIGGRAPH ’92, 1992, pp. 265–272.

[8] C. Schlick, “An inexpensive BRDF model for physically-based render-

Finally, for applications where lighting is exclusively defined ing,” in Proc. Eurographics ’94, 1994, pp. 233–246.
by environment maps, we intend to generalize our IS scheme
[9] D. Edwards, S. Boulos, P. Shirley, M. Ashikhmin, M. Stark, and
to directly process the product of the BRDF with the environ- C. Wyman, “The halfway vector disk for BRDF modeling,” ACM Trans.
ment map, following ideas presented by Jarosz et al. [40]. Graph., vol. 25, no. 1, pp. 1–18, 2006.
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 12

beige-fabric blue-metallic-paint nickel chrome

Tabulated CDFs and PDFs

Variance (0.0153, 0.0084, 0.0088) Variance (0.0006, 0.0004, 0.0033) Variance (0.0025, 0.0020, 0.0022) Variance (0.0077, 0.0062, 0.0095)
Rational Function Approximation of the inverse CDFs

Variance (0.0152, 0.0084, 0.0088) Variance (0.0009, 0.0007, 0.0048) Variance (0.0017, 0.0141, 0.0156) Variance (0.0701, 0.0631, 0.0937)
Max Lab error: 2.09 Max Lab error: 10.2 Max Lab error: 2.90 Max Lab error: 24.25
Mean Lab error: 0.43 Mean Lab error: 1.19 Mean Lab error: 0.44 Mean Lab error: 1.13
Fig. 13. Comparison of Importance Sampling with the Rational Function approximation of inverse CDFs and the original
tabulated data. 400 samples per pixel were used. The numbers under each image indicate the mean variance per color
channel. Compared to the original tabulated data where the different CDFs and PDFs require 10 MB for each material, our
RF approach requires only between 0.117 and 0.230 KB.

[10] M. Ashikhmin and S. Premoze, “Distribution-based BRDFs,” Univ. of [19] J. Kautz and M. McCool, “Interactive rendering with arbitrary BRDFs
Utah, Tech. Rep. unpublished, 2007. using separable approximations,” in Proc. EGWR ’99, 1999, pp. 247–
260.
[11] J. Wang, S. Zhao, X. Tong, J. Snyder, and B. Guo, “Modeling anisotropic
surface reflectance with example-based microfacet synthesis,” in Proc. [20] W. Matusik, H. Pfister, M. Brand, and L. McMillan, “A data-driven
ACM SIGGRAPH ’08, 2008, pp. 41:1–9. reflectance model,” in Proc. ACM SIGGRAPH ’03, 2003, pp. 759–769.
[12] K. Torrance and E. Sparrow, “Polarization, directional distribution, and
[21] J. Lawrence, S. Rusinkiewicz, and R. Ramamoorthi, “Efficient BRDF
off-specular peak in light reflected from roughened surfaces,” J. Opt.
importance sampling using a factored representation,” in Proc. ACM
Soc. Am., vol. 57, no. 9, pp. 1105–1114, 1967.
SIGGRAPH ’04, 2004, pp. 496–505.
[13] B. Cabral, N. Max, and R. Springmeyer, “Bidirectional reflection func-
tions from surface bump maps,” in Proc. ACM SIGGRAPH ’87, 1987, [22] P. Weistroffer, K. Walcott, G. Humphreys, and J. Lawrence, “Efficient
pp. 273–281. basis decomposition for scattered reflectance data,” in Proc. EGSR ’07,
2007, pp. 207–218.
[14] S. Westin, J. Arvo, and K. Torrance, “Predicting reflectance functions
from complex surfaces,” in Proc. ACM SIGGRAPH ’92, 1992, pp. 255– [23] D. Mahajan, Y.-T. Tseng, and R. Ramamoorthi, “An analysis of the in-
264. out BRDF factorization for view-dependent relighting,” in Proc. EGSR
’08, 2008, pp. 1137–1145.
[15] J. Koenderink, A. Doorn, and M. Stavridi, “BRDF expressed in terms
of surface scattering modes,” in Proc. Eur. Conf. on Comp. Vision ’96, [24] E. Lafortune, S.-C. Foo, K. Torrance, and D. Greenberg, “Non-linear
1996, pp. 28–39. approximation of reflectance functions,” in Proc. ACM SIGGRAPH ’97,
1997, pp. 117–126.
[16] P. Schröder and W. Sweldens, “Spherical wavelets: efficiently represent-
ing functions on the sphere,” in Proc. ACM SIGGRAPH ’95, 1995, pp.
[25] J. Wang, P. Ren, M. Gong, J. Snyder, and B. Guo, “All-frequency
161–172.
rendering of dynamic, spatially-varying reflectance,” in Proc. ACM
[17] T. Zickler, S. Enrique, R. Ramamoorthi, and P. Belhumeur, “Reflectance SIGGRAPH Asia ’09, 2009, pp. 133:1–10.
sharing: Image-based rendering from a sparse set of images,” in Proc.
EGSR ’05, 2005, pp. 253–265. [26] A. Ngan, F. Durand, and W. Matusik, “Experimental analysis of BRDF
models,” in Proc. EGSR ’05, 2005, pp. 117–226.
[18] M. McCool, J. Ang, and A. Ahmad, “Homomorphic factorization of
BRDFs for high-performance rendering,” in Proc. ACM SIGGRAPH ’01, [27] P. Dutré, K. Bala, and P. Bekaert, Advanced Global Illumination. A.K.
2001, pp. 171–178. Peters, 2006.
TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS, VOL. XXXX, NO. XXXXXX, XXX 20XX 13

[28] M. Kurt, L. Szirmay-Kalos, and J. Křivánek, “An anisotropic BRDF Christophe Schlick is professor in Computer Sci-
model for fitting and Monte Carlo rendering,” ACM SIGGRAPH Comp. ence at the University of Bordeaux 2 (France)
Graph., vol. 44, no. 1, pp. 3:1–15, 2010. where he has recently headed the Applied Math-
ematics and Computer Science Department. After
[29] J. Lawrence, S. Rusinkiewicz, and R. Ramamoorthi, “Adaptive numeri- having received his PhD in 1992 for his work on
cal CDF for efficient importance sampling,” in Proc. EGSR ’05, 2005, BRDF models and Monte Carlo techniques, his
pp. 11–20. research interests have embraced many aspects
[30] R. Montes, C. Urena, R. Garcia, and M. Lastra, “Generic BRDF of computer graphics, including global illumination,
sampling: A sampling method for global illumination,” in Proc. GRAPP procedural texture and geometric synthesis, curves
’08, 2008, pp. 191–198. and surfaces, point based modeling and rendering.

[31] A. Cuyt, “Recent applications of rational approximation theory: a guided

tour,” in Proc. Num. Anal. and Comp. Math. ’03, G. Psihoyios, Ed.
Weinheim, Wiley, 2003, pp. 50–52.
[32] O. Salazar Celis, A. Cuyt, and B. Verdonk, “Rational approximation of
vertical segments,” Numerical Algorithms, vol. 45, pp. 375–388, 2007.
[33] A. Wills, QPC–Quadratic Programming in C, School of Electrical
Engineering and Computer Science, University of Newcastle, Callaghan,
Australia, 2010, http://sigpromu.org/quadprog/. Xavier Granier received his Engineer and MS
[34] S. Rusinkiewicz, “A new change of variables for efficient BRDF degrees from the Grenoble Institute of Technology
representation,” in Proc. EGWR ’98, 1998, pp. 11–22. and his PhD Diploma for University Joseph Fourier
in Grenoble (France). He is currently at INRIA Bor-
[35] F. Romeiro, Y. Vasilyev, and T. Zickler, “Passive reflectometry,” in Proc. deaux Sud-Ouest (France) as a research scientist.
Eur. Conf. on Comp. Vision, 2008, pp. 859–872. His research interest includes realistic lighting and
expressive rendering, appearance modeling and
[36] A. Ben-Artzi, R. Overbeck, and R. Ramamoorthi, “Real-time BRDF acquisition. He is a member for ACM and EURO-
editing in complex lighting,” in Proc. ACM SIGGRAPH ’06, 2006, pp. GRAPHICS.
945–954.
[37] M. Stark, J. Arvo, and B. Smits, “Barycentric parameterizations for
isotropic BRDFs,” IEEE Trans. Vis. and Comp. Graph., vol. 11, no. 2,
pp. 126–138, 2005.
[38] F. X. Sillion, J. R. Arvo, S. H. Westin, and D. P. Greenberg, “A global
illumination solution for general reflectance distributions,” in Proc. ACM
SIGGRAPH ’91, 1991, pp. 187–196.
[39] R. Montes, “An importance sampling method for arbitrary BRDFs used
in global illumination applications,” Ph.D. dissertation, University of Pierre Poulin is a full professor in the Computer
Granada, 2008. Science and Operations Research department of
the Université de Montréal. He holds a Ph.D. from
[40] W. Jarosz, N. Carr, and H. Jensen, “Importance sampling spherical the University of British Columbia and a M.Sc.
harmonics,” in Proc. Eurographics ’09, 2009, pp. 577–586. from the University of Toronto, both in Computer
Science. He has served on program committees
of more than 45 international conferences. His
research interests cover a wide range of topics,
including image synthesis, image-based modeling,
procedural modeling, natural phenomena, scien-
Romain Pacanowski received his Engineer de-
tific visualization, and computer animation.
gree from EFREI and his MS degree from the
University of Bordeaux. He holds a PhD from the
University of Bordeaux (France) and Université de
Montréal (Canada). He was a postdoctoral fellow
at the CEA in 2010 and for Disney Interactive
Media Group in 2011. Since December 2011, he
is working at the CNRS as Research-Engineer.
His research interest includes realistic rendering as
well as appearance modeling and acquisition.
Annie Cuyt is full professor at the Department of
Mathematics and Computer Science of the Univer-
sity of Antwerp. She received her Doctor Scientiae
degree in 1982 from the same university, summa
cum laude and with the felicitations of the jury.
Subsequently she was a Research fellow with the
Oliver Salazar Celis received the M.Sc. degree in Alexander von Humboldt Foundation (Germany)
Computer Science in 2004 and the Doctor Scien- and she was honoured with a Masuda Research
tiae degree in 2008, both from the University of Grant (Japan). She is the author of more than 100
Antwerp, Belgium. From October 2008 until De- publications in international journals and confer-
cember 2011 he was a Post-Doctoral researcher ence proceedings, the author or editor of several
in the Department of Mathematics and Computer books and the organizer of a number of international events. Her current
Science at the same university. His research in- interests are in rational approximation theory and its applications in scientific
terests include rational approximations, numerical computing.
analysis, optimization and their use in applications
ranging from telecommunication networks to sig-
nal processing. Currently, his research focuses on
modeling and numerical aspects in computational finance.