A Data-Driven Reﬂectance Model
Wojciech Matusik ∗ Hanspeter Pﬁster † Matt Brand† Leonard McMillan ‡

Figure 1: Renditions of materials generated using our model: steel teapot with greasy ﬁngerprints (left), teapot with rust forming (right). Closeup pictures in the center. We used a spatially varying texture to interpolate between reﬂectance models for each point on the teapot.

Abstract
We present a generative model for isotropic bidirectional reﬂectance distribution functions (BRDFs) based on acquired reﬂectance data. Instead of using analytical reﬂectance models, we represent each BRDF as a dense set of measurements. This allows us to interpolate and extrapolate in the space of acquired BRDFs to create new BRDFs. We treat each acquired BRDF as a single high-dimensional vector taken from a space of all possible BRDFs. We apply both linear (subspace) and non-linear (manifold) dimensionality reduction tools in an effort to discover a lowerdimensional representation that characterizes our measurements. We let users deﬁne perceptually meaningful parametrization directions to navigate in the reduced-dimension BRDF space. On the low-dimensional manifold, movement along these directions produces novel but valid BRDFs. Keywords: Light Reﬂection Models, Photometric Measurements, Reﬂectance, BRDF, Image-based Modeling

Bidirectional Reﬂectance Distribution Functions (BRDFs) characterizes the process where light transport occurs at an idealized surface point. Traditionally, physically inspired analytic reﬂection models [Cook and Torrance 1982] [He et al. 1991] [He et al. 1992] provide the BRDFs used in computer graphics. These BRDF models are only approximations of reﬂectance of real materials. Furthermore, most analytic reﬂection models are limited to describing only particular subclasses of materials – a given model can represent only the phenomena for which it is designed. Signiﬁcant efforts have been expended on improving these models by incorporating the relevant aspects of the underlying physics. Many of these models are based on material parameters that in principle could be measured, but in practice are difﬁcult to acquire. An alternative to directly measuring model parameters is to acquire actual samples from a BRDF using some version of a goniospectro-reﬂectometer [Marschner et al. 2000] [Cornell ] [CUReT ] [STARR ] [Dana 2001] [Ward 1992] and then ﬁt the measured data to a selected analytic model using various optimization techniques [Ward 1992] [Yu et al. 1999] [Lafortune et al. 1997] [Lensch et al. 2001]. There are several shortcomings to this measure-and-ﬁt approach. First, a BRDF represented by the analytic function with the computed parameters is only an approximation of real reﬂectance; measured values of the BRDF are usually not exactly equal to the values of the analytic model. The measure-and-ﬁt approach is often justiﬁed by assuming that there is inherent noise in the measurement process and that the ﬁtting process ﬁlters out these errors. This point of view, however, ignores more signiﬁcant modeling errors due to approximations made in the analytic surface reﬂection model. Many of the salient and distinctive aspects of an objects reﬂection properties might lie within the range of these modeling errors. Second, the choice of the error function over which the optimization should be performed is not obvious. For example, error based on the Euclidean distance is a poor metric since it tends to overemphasize the importance of the specular peaks (these are usually much higher than the rest) and ignore the off-specular values. Finally, there is no guarantee that the optimization process will yield the best model. Since most BRDF models are highly non-linear, the optimization frameworks used in the ﬁtting process rely heavily on initial guesses of the models parameters. The quality of these initial guesses can have a dramatic impact on the ﬁnal

1

Introduction

A fundamental problem of computer graphics rendering is modeling how light is reﬂected from surfaces. A class of functions called
∗ MIT, Cambridge, MA. Email: wojciech@graphics.lcs.mit.edu † MERL, Cambridge, MA. Email: [pﬁster,brand]@merl.com ‡ UNC, Chapel Hill, NC. Email: mcmillan@cs.unc.edu

parameter values of the model. Another approach is to acquire dense measurements of the BRDF and use these measurements directly in the rendering process. This approach preserves those subtleties of the measured data that are lost in a data-ﬁtting approach. However, it is timeconsuming since it requires reﬂectance measurements for all desired materials in the scene. Furthermore, we end up with a collection of measured BRDFs and not with a parameterized reﬂectance model. Any change to the material property would require ﬁnding a real material with the desired property and acquiring its reﬂectance. We suggest an alternative sampling-based approach for modeling surface reﬂectance. We capitalize on the fact that it is feasible to rapidly acquire accurate reﬂectance measurements using imagebased techniques. We acquire BRDFs for a large representative set of materials. Materials in our collection include metals, paints, fabrics, minerals, synthetics, organic materials, and others. We introduce a new approach to BRDF modeling, an approach that is data driven – it interpolates/extrapolates new BRDFs from the representative BRDF data. Our approach has the advantage that the produced BRDFs look very realistic since they are based on the measured BRDFs. Furthermore, we provide a set of intuitive parameters that allow users to change the properties of the output BRDF. We also let users specify their own parameters by labeling a few representative BRDFs. We believe that this way of specifying model parameters makes our model much easier to use and control than the analytic models in which the meaning of parameters is often non-intuitive [Pellacini et al. 2000]. In our model, we do not want to store all acquired BRDFs explicitly. This leads us to the analysis of the space of all possible BRDFs for common materials in the world. A BRDF for these materials is not an arbitrary function, and we seek a representation for all possible functions corresponding to physical BRDFs. We treat each of our acquired BRDFs as a single high-dimensional vector where each measurement is an element of this vector. Then we apply both linear and non-linear dimensionality reduction tools to obtain a low dimensional manifold that characterizes the set of BRDFs we measured. In the process we also obtain a mapping between the embedding manifold and the original BRDF space. Therefore we can always compute the corresponding BRDF for each point on the manifold. An interesting side effect of our approach is that it suggests an inherent dimensionality for the space of all isotropic BRDFs. To summarize, the main contributions of this paper are: • We introduce a novel model for an isotropic BRDF that is based on measured reﬂectance for a large set of materials. • We introduce a set of perceptually-based parameters for this model. We also let users specify their own parameters. • We analyze both linear and non-linear dimensionality of the space of isotropic BRDFs. • In our model the parameter values are pre-deﬁned for many typical materials – the materials we have measured. Using our model we can also generate difﬁcult to represent effects such as rust, oxidation, or dust.

models [Cook and Torrance 1982][He et al. 1991]. An interesting transition occurred with [Ward 1992], when Ward developed a BRDF model that, while not strictly physically based, was capable of describing most signiﬁcant reﬂection phenomena. He went to great effort to ensure that his model obeyed the most basic of physical laws (reciprocity and energy conservation), and signiﬁcantly, he ﬁt his model’s parameters to actual material measurements. More recently, the availability of low-cost digital cameras has rekindled interest in BRDF acquisition and modeling. One particularly ambitious undertaking is the CUReT BRDF database [Dana et al. 1999] [CUReT ]. The CUReT database represents approximately 200 reﬂectance measurements over varying incident and reﬂected angles for a planar patch of 60 different materials. With a uniform material sample this amounts to a relatively sparsely sampled BRDF. Such a sparsely sampled BRDF is not directly useful as a table-based BRDF function; thus, it was necessary to ﬁt an analytic function in order to get a useful model. Marschner [Marschner et al. 2000] constructed another signiﬁcant BRDF measurement system. His system, although limited to only isotropic BRDF measurements, was both fast and robust. In particular, his system took unique advantage of reciprocity, bilateral symmetry, and multiple simultaneous measurements to achieve unprecedented leverage from each reﬂection measurement. This offers a signiﬁcant advantage. It ﬁlters measurement noise due to minute variations over the surface, errors due to spatial variations in photosite response within the image sensor, and variations in illumination intensity. In the face of such statistical averaging, one is hard pressed to attribute the inevitable residual errors that occur when model ﬁtting to additional systematic noise, rather than failings of the analytic model. The inherent dimensionality of a BRDF, combined with the desire to sample it at high resolutions in order to model specular, incident, and retroreﬂection effects, leads to an unwieldy sampling and storage problem. Many researchers have addressed this speciﬁc problem by searching for a more appropriate basis for representing BRDFs. Spherical harmonics [Westin et al. 1992] and spherical wavelets [Schr¨ der and Sweldens 1995] are natural choices for o representing the angular parameters of the BRDF. Other efﬁcient representations include wavelets [Lalonde and Fournier 1997], Zernicke polynomials [Koenderink et al. 1996],and separable approximations obtained using singular value decomposition [Kautz and McCool 1999] or a purely positive matrix factorization [McCool et al. 2001]. Furthermore, recent image-based approaches to BRDF modeling [Lensch et al. 2001] have demonstrated the power of using linear combinations of a compact reﬂectance function basis sets for modeling spatially varying BDRFs. Such linear decompositions lead to an interesting question: can the true space of potential BRDFs be described as a linear combination of basis functions? Clearly, factorizations of the sort used to compress BRDFs are linear, allowing for arbitrary mixtures of their basis vectors to ﬁt a given set of data. If this decomposition approach were in fact valid, it would imply that linear combinations of actual BRDFs might be used to model original and physically plausible reﬂection models. Exploring the ramiﬁcations of this hypothesis is one of our motivations for developing a sample-based generative model.

2

Previous Work

3

Data Acquisition

The value of physically accurate reﬂectance models has long been understood within the computer graphics community [Blinn 1977]. The availability of BRDF models based on the actual physics of light transport and validated by empirical measurements were a signiﬁcant catalyst in this realization [Torrance et al. 1966] [Trowbridge and Reitz 1975]. Physical accuracy was an impetus behind the development of many subsequent computer graphics reﬂection

In order to acquire a sufﬁcient number of adequately sampled BRDFs, it was necessary to build a measurement device. Our modeling approach placed two requirements on the acquired data: ﬁrst, that each BRDF be sampled densely enough that it could be used directly as a table-based model, and second, that the space of BRDFs be sampled adequately so as to span the range of models that we hope to generate. Accordingly, we have built a BRDF measurement device suitable for rapidly acquiring high-quality BRDFs for a wide

is not possible to reproduce original images (the specular highlight becomes an oval shape, oriented at different directions). We use a different coordinate system, described in [Rusinkiewicz 1998] and shown in Figure 3. This coordinate frame is based on the angles with respect to the half-angle (half-vector between incoming and outgoing directions). This coordinate frame allows us to vary the sampling density near the specular highlight. Speciﬁcally, we vary θh (angle between the normal and the half-vector), assigning smaller bins for values near specular reﬂection and larger bins for angles far away from the specular reﬂection.

Figure 2: A photograph of our high-speed BRDF measurement gantry. range of different materials (see Figure 2 ). The image-based BRDF measurement device described by [Marschner et al. 2000] inspired our design. Our acquisition system requires a spherically homogenous sample of the material. The system is placed in a completely isolated room painted in black matte. It consists of the following components: a QImaging Retiga 1300 (a 10 bit, and a 1300x1030 resolution Firewire camera), a Kaidan MDT-19 (a precise computercontrolled turntable), and a Hamamatsu SQ Xenon lamp (a lamp with stable light emission output and a continuous and relatively constant radiation spectrum over the visible light range). The lamp is mounted on an arm to the turntable. The light orbits the measurement sample placed at the center of rotation; the camera is stationary. Our camera is geometrically calibrated using the technique described in [Zhang 1998]. The position of the light source is determined using a contact digitizer (FARO Arm). We use the same digitizer to determine the position of the center of the material sample. The radius of the sample is measured with calipers. The light source moves in increments of approximately 0.5 ◦ from the point exactly opposite the camera (the sample is in between the camera and the light source) to the point exactly in front of the camera. We take a total of 330 high dynamic range pictures to cover the required half circle. This process takes about 3 hours. For each high dynamic range picture we take a total of 18 10-bit photographs. The exposure time ranges from 40 microseconds to 20 seconds. We use the fact that our CCD camera has a very linear response curve to derive the high dynamic measurement. For each pixel in the image we ﬁt a line to the exposure time vs. radiance values. The slope of the line is used as the radiance estimate. The correlation of this line is higher than 0.998. Each acquired image of the sample sphere represents many BRDF samples. Essentially each pixel of the sphere is treated as a separate BRDF measurement. In order to compute the speciﬁc BRDF value for a given pixel we perform the following steps. First, we intersect the ray deﬁned by the pixel with the sphere to determine point P. Then, we compute the normal at point P on the sphere, the vector and the distance to the light source, and the vector to the camera pixel. Next, we compute the irradiance at point P due the light source (taking into account distance to the light source and foreshortening). Finally, we compute the BRDF value as the ratio of the high dynamic range radiance to the irradiance.

Figure 3: The standard coordinate frame is shown on the left. Rusinkiewicz’s coordinate system is shown on the right. We still discretize θh ,θd into 90 bins and φd into 360 bins. This results in a total of 90 x 90 x 360 = 2,916,000 bins for each color component. We halve this number to 1,458,000 by enforcing the reciprocity constraint: f (θh , θd , φd ) = f (θh , θd , φd + π ) With this constraint we need only to discretize φd into 180 bins. (1)

Figure 4: Two log images of a sphere (alumina oxide). A real image is shown on the left. A synthesized image using tabulated BRDF data is shown on the right. Our measurement process gives us typically 20-80 million BRDF samples for each material. We reduce the noise in our measurements by removing the outliers in each bin (lowest and highest 25% of the values), and we average the remaining measurements. This statistical smoothing is intended to remove systematic noise as well as compensate for small variations in material properties over the sample. As a ﬁnal validation we render a synthesized version of our sample sphere and compare it to the corresponding acquired high dynamic range image. We conduct this inspection for all input light conﬁgurations. Pictures for a typical acquired material are in shown in Figure 4. The rendered images reproduce the input images very well. We have used our device to acquire BRDF measurements of more than 130 different materials, including metals, plastics, painted surfaces, and cloth. Figure 5 depicts some of the materials that were sampled. We have removed from further analysis some materials that exhibited signiﬁcant subsurface scattering, anisotropy, or non-homogenity.

4

Data Representation

We found that specular peaks were difﬁcult to represent using the natural coordinate system (θin ,θout , φdi f f ). Even when binning a BRDF at a dense grid (every 1 ◦ spacing for each dimension), it

Figure 6: Rendered teapots using BRDFs from our database: nickel, hematite, gold paint, and pink fabric.

Figure 5: Pictures of 100 of our acquired materials.

5

Data Analysis

These sampled BRDFs can be used directly by a renderer. Several examples of that are shown in Figure 6, where a teapot is rendered under natural illumination using the raw acquired data. Our ultimate goal, however, is to construct an empirical BRDF model that can be used to generate novel, yet plausible, reﬂectance functions directly from this database. We begin with the following assumption: if we treat each of our BRDF samples as a high dimensional vector in an abstract BRDF space, we expect that all physical BRDFs lie upon a lower dimensional manifold within this space indicative of their inherent dimensionality. This is a common assumption used by others [Cula and Dana 2001] and it is consistent with the relatively small number of parameters seen in analytic BRDF models. Therefore, we breakdown the task of constructing an empirical BRDF model into two phases: discovering this lower dimensional model, and deﬁning an interpolation scheme within this lower-dimensional subspace.

Figure 7: Plot of the eigenvalues resulting from PCA of the data set. We perform the analysis in the log space (we apply the natural logarithm to each element of vector X). There are several reasons for this normalization. First, there is a huge difference (on the order of a few magnitudes) between the specular and non-specular values of the BRDF. If used in the original space, the analysis tools would associate more importance to noise in the specular values than the actual non-specular components. The linear analysis would depreciate importance of these non-specular values (the non-specular values are perceptually important). Our operation is also justiﬁed by the fact that the human visual system is sensitive to ratios rather than absolute radiance values. Singular value decomposition was then applied to X T X (a 104x104 matrix). The singular values in this case are the squares of the desired eigenvalue magnitudes. A plot of these eigenvalues is shown in Figure 7. We also show in Figure 8 the reconstruction of a typical material using ﬁrst 1, 5, 10, 20, 30, 45, 60, and all principal components. We see that good reconstruction is usually obtained using the ﬁrst 30-40 components. While there is a considerable fall off in the sequential values seen in this plot, the plateau is reached around 45th eigenvalue (the reconstruction error is about 1% at that point). This dimension of the embedding subspace is considerably higher than our intuition would suggest, based on the typical number of parameters used in analytic BRDF models. We veriﬁed that the 45-dimensional space deﬁned by the ﬁrst principal components reconstructs all our measured BRDFs well. However, it spans a space that is bigger than the space of all possible BRDFs. We are able to ﬁnd the points in this

5.1

Linear Analysis

In the case where the physical BRDF manifold lies on a linear subspace, the analysis tools for both manifold discovery and interpolation are well known. In this case, Principal Component Analysis (PCA)[Bishop 1995] effectively determines a set of basis vectors that span the desired subspace, and linear combinations of samples can be used for interpolation. Linear manifold approaches have proven extremely effective in some problem domains, such as face synthesis [Blanz and Vetter 1999] and radiance interpolation [Chen et al. 2002]. Potential linear manifolds are generally suggested when there is a noticeable plateau in the magnitudes of the sorted eigenvalues. When this plateau occurs on the kth eigenvalue, we can model the data as a k-dimensional linear subspace with a residual error bounded by the square root of the sum of the squares of the remaining eigenvalues. We began our analysis of the BRDF samples by searching for a linear embedding manifold (a hyperplane). The three color channels of each BRDF sample were assembled into a column vector and concatenated to form a 4,374,000 by 104 measurement vector matrix X.

Figure 8: Reconstruction of a BRDF from principal components in the order of increasing number of components – mean, 5, 10, 20, 30, 45, 60, and all. subspace that do not correspond to any physical materials. In other words, using linear combinations of the components, we can obtain the data samples that do not look like BRDFs. We illustrate this point in Figure 9. Moreover, in order to span the whole space, we would need to have at least 45 parameterization directions in order to reach all speciﬁed BRDFs. This suggests that the space of all possible BRDFs lies on a lower-dimensional manifold that is nonlinearly embedded in the 45D linear space. In the next section, we apply recently developed nonlinear dimensionality reduction techniques to discover this lower dimensional manifold.

5.2

Nonlinear dimensionality reduction
Figure 10: A simple charting example. Points (⊕) sampled from a unknown manifold (gray curve) are projected onto three subspaces (red, green, and blue lines) and assigned a probability (indicated by size) according to their distance from the point where the chart touches the manifold. A minimal-distortion merger of these charts gives a ﬂattening of the manifold in a lower dimensional space, where the mapped locations of points are the probability-weighted combinations of their chart-speciﬁc locations. cover the data manifold, in the sense that adjoining Gaussians have similar orientation. The dominant axes of each Gaussian specify a subspace. Projecting the data into this subspace gives a “chart” of one part of the manifold. A chart preserves local structure where it touches the manifold and suppresses measurement noise that displaces samples off the manifold. A data point has a location and a probability in every chart. Due to curvature of the manifold, a chart gives a very distorted picture of faraway points; these points are assigned very low probability. The pancake Gaussians are solved under a criterion that optimizes the charts for the ensuing “connection.” The connection is an afﬁne merger of all charts in the target space–effectively a ﬂattening of the manifold that minimally distorts all charts and maximizes agreement between overlapping charts of the locations of points to which they assign high probability. The connection gives

Nonlinear dimensionality reducers (NLDR) compute low-distortion embeddings of high-dimensional data in low-dimensional target spaces. The nonlinearity usually obtains from the fact that only local relationships in the ambient space are preserved while longdistance relationships are presumed to be corrupted by the curvature of the manifold in the ambient space. First-generation NLDRs such as nonmetric MDS [Kruskal and Wish 1978], IsoMap [Tenenbaum et al. 2000], and LLE [Roweis and Saul 2000] generalize PCA to give low-dimensional embeddings of the data, but offer no mapping of the data points. Recently, two second-generation methods have been announced that offer continuous mappings between the embedding an the original (ambient) space: Automatic Alignment [Teh and Roweis 2003] combines LLE with a set of pre-estimated local dimensionality reducers–each of which is presumed to be ﬁtted to a relatively ﬂat subset of the manifold–and solves for a mixture of these projections that globally ﬂattens the data while minimizing barycentric distortion in each point neighborhood. Charting [Brand 2003] solves for a kernel-based mixture of projections that minimizes Euclidean distortion of local neighborhoods; it includes a solution for the local dimensionality reducers needed by automatic alignment. We chose to use charting because it is explicitly designed to work well with small numbers of samples and to suppress measurement noise, two conditions that tend to break methods for dimensionality reduction from local relationships. Figure 10 gives the main geometric intuition behind charting. First one solves for a set of ﬂat “pancake” Gaussians that smoothly

Figure 9: Nonlinear spaces generate valid BRDFs where linear spaces fail. Original BRDF corresponding to a point A on a 45 dimensional hyperplane (left). Physically implausible reﬂectance (hole in the middle of the specular highlight) corresponding to moving away from a point A on the 45 dimensional linear subspace (center). Physically plausible reﬂectance corresponding to moving equally far away from point A on the 15 dimensional non-linear manifold (right). would ﬂatten the cone into an annulus. Each ﬂattening would suppress the noise in different directions, some more fortuitous than others. While the 10D manifold exhibits good reconstruction of the original data, our goal is to synthesize novel BRDFs. With that in mind we chose to work on a 15D manifold because interpolations on it pass even closer to the data density (with error comparable to 45D PCA reconstruction). Moreover, this dimensionality is roughly consistent with previous isotropic BRDF models [Ward 1992], [Lafortune et al. 1997], and [Koenderink et al. 1996], which have at least 10 degrees of freedom. A charted manifold of BRDF data makes it possible to treat the space of BRDFs as if it were linear, and to identify meaningful axes of variation in this embedding space. An interpolating or extrapolating line in this space is a nonlinear curve in the original BRDF space that passes closer to the data density that the equivalent straight line would (on average), simply because it stays on the manifold where a straight line does not. This translates directly to superior BRDF synthesis, as will be demonstrated below. Figure 11: Data reconstruction error as a function of the dimensionality of the global chart. The sharp drop in this error curve indicates that a 10-dimensional chart is sufﬁcient for the BRDF data. In fact, that chart has a better reconstruction error than a 25dimensional PCA. mappings between the ambient and target spaces, which are simply mixtures of afﬁne projections, weighted by the probability that a point “belongs” to each chart. The dimensionality-reducing mapping from the ambient to the target space effectively imposes a lowdimensional coordinate system on the samples, while the inverse mapping gives a smoothly curving low-dimensional surface in the ambient space, effectively reconstructing the original manifold. For charting, one must specify a set of chart centers, a width parameter σ for the Gaussians, and a target dimensionality d. We used the default settings: one chart centered on each data point and σ = half the average distance between each point and its closest neighbor. Note that locating a chart on a point does not cause the manifold to pass through that point–only near it. See [Brand 2003] for additional details. As with PCA, the data-reconstruction error of a charted data set gives an indication of the true dimensionality of the manifold. Figure 11 shows that our BRDF data probably lies on a 10D manifold. The reconstruction error does not decline monotonically because each dimensionality may merit a different ﬂattening. For example, if the data were sampled from a truncated cone, the best 1D chart would simply be height along the cone, while the best 2D chart

6

Model Construction

In order to use our sample-based reﬂectance model it is necessary to develop intuitive user interfaces for specifying and exploring new materials. We investigated methods for characterizing material traits by analogies derived from the existing samples. We believe that such methods provide the best and most intuitive user interface [Pellacini et al. 2000]. Our model is built from actual physical measurements and it reproduces these measurements. Therefore, we have deﬁned model parameters for a large collection of materials – materials we have measured. We believe that the most useful scheme of navigation is when users can choose as a starting point some type of the material similar to the one they desire. In our case they can pick any of the measured materials. Then, they would change the reﬂectance properties of this material according to one of the following schemes (these navigation schemes are applicable for both linear and nonlinear manifold models). The simplest method is to choose another BRDF and move in this direction. Although of limited use, this method works well for perceptually similar materials. A more useful approach is to deﬁne directions corresponding to a desired trait (the parameterization direction is a 45D vector for linear space and a 15D vector for nonlinear space). We pick some arbitrary point on the manifold and then move in the direction deﬁned by the vector by adding it to the current position to increase the trait, or subtracting it to decrease the trait. We can backproject the current point onto the original BRDF space to check the corresponding BRDF. Next, we describe various procedures for identi-

fying trait vectors. Our modeling approach requires the user to specify a sufﬁcient set of traits. This speciﬁcation can be as simple as a binary classiﬁcation (i.e., noting whether each acquired BRDF has the speciﬁed trait.) We also allow the user to leave a BRDF unspeciﬁed in cases where the trait is hard to determine or simply does not apply. Usually the more samples we specify for each class the more precise the direction is. There are many different ways to deﬁne the parameterization directions based on the classiﬁcation. We have examined and evaluated a few. (A) Mean difference [Blanz and Vetter 1999]: In this approach we compute the average of each BRDF in each complementary pair of clusters associated with a trait (i.e., those samples with, and those without) in the embedding space. Then the vector between these complement averages in the embedding space is the parameterization direction. This direction vector is then applied (added or subtracted) to the current point in the embedding space. (B) Support vector machines [Vapnik 1995]: Support vector machines determine the hyperplane which separates the data points in the ﬁrst material class from the data points in the second class. The partitioning hyperplane has maximum distance to the closest points (called support vectors) in both material classes. The parameterization direction is the normal to this hyperplane. The hyperplane is deﬁned in 15D space for non-linear analysis and 45D for the linear space. This method also tells us on which side of the hyperplane the current point is, and how far the point is from the plane. (C) Fisher’s linear discriminant [Duda and Hart 1973]: Each material class corresponds to a some distribution of high-dimensional data (15D for non-linear analysis and 45D for linear analysis). Fisher’s linear discriminant deﬁnes a projection of these distributions on the axis such that the distributions projected on this axis are the most separable (the projection maximizes the distance between the means of the two classes while minimizing the variance of each class). In practice, support vector machines performed the best on our data set and Fisher’s linear discriminant performed the poorest. Since we want our model to preserve the basic principles of physics, we have to disallow movements on the manifold that do not adhere to these principles. We consider the three following principles: • Reciprocity: As mentioned before, reciprocity in our model is met by default since we store only half of the BRDF vector. • Non-negativity: We allow the user to move only in the space so that all the values in the backprojected vector are positive. • Energy conservation: A unit of light energy is applied at some incoming light direction. If the sum of energy in all outgoing directions is less than one (we assume that the surface does not emit energy by itself) then the energy is conserved. This has to be true for all incoming light directions in order for a BRDF to follow energy conservation. We enforce this and do not allow the users to produce BRDFs for which the sum of energy for any incoming direction is greater than one.

Figure 12: Diffuseness trait vs specularness trait. Observe that the diffuseness and specularness traits exhibit a weak inverse correlation. The green, blue, and red vectors denote projections of the BRDF interpolations shown in the second, third, and fourth rows of Figure 16 respectively.

Figure 13: Metallic-like trait vs specularness trait. Observe that the metallic-like and specularness traits exhibit a weak correlation. The green, blue, and red vectors denote projections of the BRDF interpolations shown in the second, third, and fourth rows of Figure 16 respectively. metallic-like, plastic-like, roughness, silverness, gold-like, fabriclike, acrylic-like, greasiness, dustiness, rubber-like. In a sense, these parameters are arbitrary since the classiﬁcation is completely based on the subject’s interpretation. We could have chosen traits without physical connotations, such as ugly or pleasing. Alternatively, the traits could have been based on actual measurable quantities, such as conductivity and mean surface variation. Our test subject characterized each BRDF as one of three choices: 1) possessing the particular trait, 2) not possessing the trait, or 3) unclear. We then used the subject’s characterizations to build a model in both the linear and non-linear embedding spaces using Support Vector Machines. The results from this trait-based analysis are shown as projections onto the derived trait vectors in Figures 12, 13, and 14. These projections are computed in the linear embedding space given by

7

Results

Once the BRDFs are acquired and validated, as described in section 4, we performed both linear and non-linear dimensionality reduction as described section 5. We then set out to construct a perceptual BRDF model using the techniques outlined in section 6. This section presents the results from a typical model construction session. A test subject was asked to characterize each of the BRDFs from our database using 16 different traits. These included redness, greenness, blueness, specularness, diffuseness, glossiness,

our non-linear model. Observe that the metallic and specular characteristics are weakly correlated, the specular and diffuse traits are weakly inverse-correlated, and the glossy and diffuse traits are inverse-correlated. This is what we would expect. Note that we do not make attempts to model independent traits in either our trait selection or trait vector derivations. Therefore, we expect that addition of a particular trait to an existing BRDF may effect other traits. This lack of parameter independence is a tradeoff that we accept in order to establish perceptually meaningful parameters in our modeling approach. Despite the fact that the parameterization vectors are not orthogonal, they did span the whole 15D non-linear embedding space and provide an intuitive set of “dials” for users to design materials.

metallic trait is correlated to the glossiness and specularness traits used as axis. The fourth row starts with YellowDiffusePaint BRDF and shows the addition of the glossiness trait, which is depicted as the red path in Figures 12, 13, and 14. The direction of this path is as we would expect, and it has a large magnitude due to the fact the the YellowDiffusePaint BRDF is located far away from the glossy examples in the projections shown. Overall the non-linear basis set results in a more robust model than our linear basis set, in that we were able to move large distances within the non-linear embedding space and still generate physically plausible BRDFs with the expected appearance. Our modeling approach also allows us to associate approximate trait vectors with physical processes. This can be done in one of two ways, by ﬁtting a least-squares line to a path of speciﬁed BRDFs in the embedding space, or by computing local piecewise difference vectors between examples. As an example, we have modeled metal oxidation. We measured the reﬂectance changes as a metal was exposed to an acidic environment. It changed from highly specular polished material to black matte material. The acquired four BRDFs determine a path in the embedding space. The intermediate stages are interpolated in the embedding space and backprojected to the sample space (Figure 17). Figure 1 illustrates another process – rust formation. We used a spatially varying texture to select rust levels for each point on the teapot. We are currently measuring more processes like this such as copper patination and other types of rust formation.

8

Future Work and Conclusions

Figure 14: Glossiness trait vs diffuseness trait. Observe that the glossiness and diffuseness traits exhibit an inverse correlation. The green, blue, and red vectors denote projections of the BRDF interpolations shown in the second, third, and fourth rows of Figure 16 respectively. Once trait vectors are established, we can add and subtract them from our data points in our embedding space. In Figure 15 we demonstrate four examples of varying user-speciﬁed traits using the linear model. The ﬁrst row shows a teapot rendered using our GoldPaint BRDF on the far left, and the effect of adding the redness trait in successive steps to the right. The second row starts from a our SpecularGold BRDF (left) with successive additions of the silverness trait. The third row adds the gold-like trait to the BlueGlossyPaint BRDF (left). Finally, the fourth row shows the addition of the specularness trait to the BlackMattePlastic BRDF. It is our experience that the linear model gives reasonable BRDFs if only small displacements are permitted. If the displacement is too large, physically invalid BRDFs result (as illustrated in Figure 9). We then applied the same trait classiﬁcations and Support Vector Machine calculations to the embedding space of our non-linear model. Figure 16 demonstrates 4 examples using this approach. The ﬁrst row of Figure 16 shows our Copper BRDF on the left, with successive additions of the roughness trait. The second row begins with our GreenAcrylic BRDF and shows the addition of the blueness trait. The trajectory of this path is also illustrated in Figures 12, 13, and 14. Notice that color-change speciﬁcation is not particularly correlated with any of traits used for these projections. Thus, we would expect relatively small movements and no preferred direction. The third row, on the other hand, represents the addition of the metallic trait to the VioletAcryllic BRDF model, whose path is also illustrated in the projections in Figures 12, 13, and 14. The path trajectory of this example conforms to our expectations, and its magnitude is large in these visualizations since the

In this paper we have introduced a new approach for modeling isotropic BRDFs. Our model generates new surface reﬂectance models by forming combinations from a set of densely-sampled, acquired BRDFs. We are hopeful that data-driven reﬂectance modeling approaches, like ours, can greatly expand the range of material models used in computer graphics rendering. In order to develop an effective and efﬁcient interpolation scheme we choose to ﬁrst analyze the inherent dimensionality of our data set. To this end we applied both in linear subspace and nonlinear manifold analysis. The results of this analysis are suggestive of the overall structure of BRDFs. Speciﬁcally, we found that the linear subspace model lent itself to the creation of physically implausible BRDFs, and a large number of dimensions (around 45) were required to adequately represent our measurements. Nonetheless, we still found the linear subspace model to be useful for interpolation over small distances. The nonlinear model, on the other hand, was much more compact in its dimensionality (around 14 dimensions for the same accuracy as the 45-dimension linear subspace model), and more robust in its ability to interpolate plausible BRDFs over long distances. However, we caution against over generalizing from our results. We are comfortable in saying that our modeling approach effectively represents our data set, but our sample size is still relatively small to draw conclusions regarding the fundamental nature of isotropic BRDFs. However, we are optimistic that techniques like ours can be used to greatly expand our knowledge in these areas. We also have demonstrated methods for deﬁning intuitive parameters for navigating within BRDF models. These techniques can easily be customized for a range of industrial and artistic applications. Furthermore, they can be personalized for individual use or made objective by incorporating physical measurements. The advantages of our data-driven BRDF model include a high degree of realism, a perceptually meaningful parameterization, relative ease of modeling for complex surface materials, and speed of evaluation. The main disadvantage of the model is its size. We believe that the model we propose can easily be incorporated into existing rendering systems. We also hope to extend our

Figure 15: Navigation in the linear space. Each row corresponds to changing one parameter of the model. The ﬁrst row shows an increase in the redness trait applied to the GoldPaint BRDF. The second row illustrates an increase in the silverness trait applied to the SpecularGold BRDF. Row three applies the gold-like trait to the BlueGlossyPaint BRDF. The fourth row shows an increase in the specularness trait applied to the BlackMattePlastic BRDF. work in sample-based reﬂectance modeling to include anisotropy (4D BRDF), macro-scale surface variations typically described by BTFs, and subsurface scattering effects (BSSRDF). Another obvious extension would be to use this model in solving inverse rendering problems.
B RAND , M. 2003. Charting a manifold. In To appear, Proc. NIPS-15. C HEN , W.-C., B OUGUET, J.-Y., C HU , M. H., AND G RZESZCZUK , R. 2002. Light Field Mapping: Efﬁcient Representation and Hardware Rendering of Surface Light Fields. ACM Transactions on Graphics 21, 3 (July), 447–456. ISSN 0730-0301 (Proceedings of ACM SIGGRAPH 2002). C OOK , R., AND T ORRANCE , K. 1982. A reﬂection model for computer graphics. ACM Transactions On Graphics 1, 1 (January), 7–24. C ORNELL. Light Measurement Laboratory. http://www.graphics.cornell.edu/ research/measure/. Web page.

9

Acknowledgments

We would thank Steven Gortler and Julie Dorsey for helpful discussions; Fredo Durand for help with writing the shader; Henrik Wann Jensen for the Dali renderer; Joe Marks for his continuing support; the anonymous reviewers for their constructive comments; and Jennifer Roderick Pﬁster for proofreading the paper. Parts of this work were supported by NSF CAREER grant 9875859.

C ULA , O. G., AND DANA , K. 2001. Compact representation of bidirectional texture functions. In Proc. of Computer Vision and Pattern Recognition, in press. CUR E T. Columbia-Utrech Reﬂectance and Texture. http://www.cs.columbia.edu/ CAVE/curet/. Web page.

References
B ISHOP, C. M. 1995. Neural Networks for Pattern Recognition. Clarendon Press. B LANZ , V., AND V ETTER , T. 1999. A morphable model for the synthesis of 3D faces. Computer Graphics 33, Annual Conference Series, 187– 194. B LINN , J. 1977. Models of light reﬂection for computer synthesized pictures. Computer Graphics 11, Annual Conference Series, 192–198.

DANA , K. J., VAN G INNEKEN , B., NAYAR , S. K., AND KOENDERINK , J. J. 1999. Reﬂectance and texture of real world surfaces. ACM Transactions on Graphics 1, 18, 1–34. DANA , K. 2001. BRDF/BTF measurement device. ICCV, 460–466. D UDA , R. O., AND H ART, P. E. 1973. Pattern Classiﬁcation and Scene Analysis. Wiley-Interscience Publication, New York. H E , X., T ORRANCE , K., S ILLON , F., AND G REENBERG , D. 1991. A comprehensive physical model for light reﬂection. Computer Graphics 25, Annual Conference Series, 175–186.

Figure 16: Navigation on the non-linear manifold. Each row corresponds to changing one parameter of the model. The ﬁrst row shows an increase in the roughness trait applied to the Copper BRDF. The second row illustrates an increase in the blueness trait applied to the GreenAcrylic BRDF. Row three applies the metallic-like trait to the VioletAcryllic BRDF. The fourth and ﬁfth rows show an increase in the glossiness trait applied to the YellowDiffusePaint BRDF (the images in the ﬁfth row are the enlarged versions of the ﬁrst and last image in the fourth row).

Figure 17: Progression of the steel oxidation process. Our BRDF model also supports interpolation along physically meaningful path. In this example we start with a completely polished steel sample (upper left) and progressively oxidize it. The ﬁnal black-oxidized sample is shown on the lower right.
H E , X., H EYNEN , P., P HILLIPS , R., T ORRANCE , K., S ALESIN , D., AND G REENBERG , D. 1992. A fast and accurate light reﬂection model. Computer Graphics 26, Annual Conference Series, 253–254. K AUTZ , J., AND M C C OOL , M. 1999. Interactive rendering with arbitrary BRDFs using separable approximations. In Rendering Techniques ’99 (Proceedings of the Tenth Eurographics Workshop on Rendering), Springer Wien, New York, NY, 281–292. KOENDERINK , J., VAN D OORN , A., AND S TAVRIDI , M. 1996. Bidirectional Reﬂection Distribution Function expressed in terms of surface scattering modes. European Conference on Computer Vision. K RUSKAL , J. B., AND W ISH , M. 1978. Multidimensional Scaling. Sage Publications, Beverly Hills, CA. L AFORTUNE , E., F OO , S.-C., T ORRANCE , K., AND G REENBERG , D. 1997. Non-linear approximation of reﬂectance functions. Computer Graphics 31, Annual Conference Series, 117–126. L ALONDE , P., AND F OURNIER , A. 1997. A wavelet representation of reﬂectance functions. IEEE Transactions on Visualization and Computer Graphics 3, 4, 329–336. L ENSCH , H., K AUTZ , J., G OESELE , M., H EIDRICH , W., AND S EIDEL , H.-P. 2001. Image-based reconstruction of spatially varying materials. In Proceedings of the 12th Eurographics Workshop on Rendering, 104– 115. M ARSCHNER , S., W ESTIN , S., L AFORTUNE , E., AND T ORRANCE , K. 2000. Image-based measurement of the Bidirectional Reﬂection Distribution Function. Applied Optics 39, 16. M C C OOL , M., A NG , J., AND A HMAD , A. 2001. Homomorphic factorization of BRDFs for high-performance rendering. Computer Graphics 35, Annual Conference Series, 171–178. P ELLACINI , F., F ERWERDA , J., AND G REENBERG , D. 2000. Toward a psychophysically-based light reﬂection model for image synthesis. Computer Graphics 34, Annual Conference Series, 55–64. ROWEIS , S. T., AND S AUL , L. K. 2000. Nonlinear dimensionality reduction by locally linear embedding. Science 290 (December 22), 2323– 2326. RUSINKIEWICZ , S. 1998. A new change of variables for efﬁcient BRDF representation. In Rendering Techniques ’98 (Proceedings of Eurographics Rendering Workshop ’98), Springer Wien, New York, NY, G. Drettakis and N. Max, Eds., 11–22. ¨ S CHR ODER , P., AND S WELDENS , W. 1995. Spherical wavelets: Efﬁciently representing functions on the sphere. Computer Graphics 29, Annual Conference Series, 161–172. STARR. NIST reference reﬂectometer: STARR facility. Web page. http://physics.nist.gov/Divisions/Div844/facilities/brdf/starr.htm. T EH , Y. W., AND ROWEIS , S. T. 2003. Automatic alignment of hidden representations. In To appear, Proc. NIPS-15. T ENENBAUM , J. B., DE S ILVA , V., AND L ANGFORD , J. C. 2000. A global geometric framework for nonlinear dimensionality reduction. Science 290 (December 22), 2319–2323. T ORRANCE , K., S PARROW, E., AND B IRKEBAK , R. 1966. Polarization, directional distribution, and off-specular peak phenomena in light reﬂected from roughened surfaces. Journal of the Optical Society of America 56, 7 (July), 916–925. T ROWBRIDGE , T., AND R EITZ , K. 1975. Average irregularity representation of roughened surfaces. Journal of the Optical Society of America 65, 5, 531–536. VAPNIK , V. 1995. The Nature of Statistical Learning Theory. Springer, N.Y. WARD , G. 1992. Measuring and modeling anisotropic reﬂection. Computer Graphics 26, Annual Conference Series, 265–273. W ESTIN , S., A RVO , J., AND T ORRANCE , K. 1992. Predicting reﬂectance functions from complex surfaces. Computer Graphics 26, Annual Conference Series, 255–264. Y U , Y., D EBEVEC , P., M ALIK , J., AND H AWKINS , T. 1999. Inverse global illumination: Recovering reﬂectance models of real scenes from photographs. In Computer Graphics, SIGGRAPH 99 Proceedings, 215– 224. Z HANG , Z. 1998. A ﬂexible new technique for camera calibration. Tr98-71, MSR Redmond.