Compact and intuitive data-driven BRDF models

Measured materials are rapidly becoming a core component in the photo-realistic image synthesis pipeline. The reason is that data-driven models can easily capture the underlying, fine details that represent the visual appearance of materials, which can be difficult or even impossible to model by hand. There are, however, a number of key challenges that need to be solved in order to enable efficient capture, representation and interaction with real materials. This paper presents two new data-driven BRDF models specifically designed for 1D separability. The proposed 3D and 2D BRDF representations can be factored into three or two 1D factors, respectively, while accurately representing the underlying BRDF data with only small approximation error. We evaluate the models using different parameterizations with different characteristics and show that both the BRDF data itself and the resulting renderings yield more accurate results in terms of both numerical errors and visual results compared to previous approaches. To demonstrate the benefit of the proposed factored models, we present a new Monte Carlo importance sampling scheme and give examples of how they can be used for efficient BRDF capture and intuitive editing of measured materials.


Introduction
Over the last decade, we have seen the development of computer graphics algorithms and techniques which enable the quality and accuracy needed to make rendered images truly comparable to photographs of the same scene.The accuracy in the simulation of scattering at surfaces and computation of the light transport in a scene is determined by the way the material properties such as color and reflectance, modeled by the Bidirectional Reflectance Distribution Function (BRDF) [30], are measured and represented.This has led to the development of both accurate parametric models [3,8,11,16,24,35,47], and advanced data-driven methods [20,26,29] for describing, measuring and analyzing for virtually all classes of materials.
B Tanaboon Tongbuasirilai tanaboon.tongbuasirilai@liu.se Jonas Unger jonas.unger@liu.seMurat Kurt murat.kurt@ege.edu.tr 1 Linköping University, Linköping, Sweden 2 International Computer Institute, Ege University, Izmir, Turkey Recently, measured materials and appearance capture techniques have gained significant popularity, for an overview see [14,23].The advantage of data-driven approaches is that carefully measured BRDFs automatically bring the detailed appearance of real-world materials into rendering pipeline.There are, however, still a number of important challenges that need to be solved in order to make data-driven models flexible and easy to use in practice.First, although significant progress has been made, [1,31], accurate BRDF and SvBRDF measurements are still challenging due to mechanical and computational complexity.Second, measured data are difficult to edit, analyze, and in other ways interact with.It is, therefore, necessary to further develop compact representations which allow for detailed representations of the reflectance distributions with minimal approximation errors.
In this paper, we set out to address these challenges by developing a new class of separable data-driven BRDF models which allow for compact representations of BRDFs in a low-dimensional space.We focus on isotropic BRDFs and propose a new separable BRDF representation inspired by the projected deviation vector (PDV) parameterization described by Löw et al. [24].The result is a compact and intuitive data-driven model, in which BRDFs can be described using independent one-dimensional (1D) strictly nonneg-ative factors, with a smaller error than previous factored representations.The independent 1D factors also enable efficient importance sampling, which is important to Monte Carlo rendering algorithms.In contrast to previous BRDF factorization methods, see for example [7,38], our separable BRDF model is described by intuitive factors which enable detailed user control for analysis or editing.The main contributions of this paper are: -A study of the separability of existing data-driven BRDF models.-Two new BRDF models specifically designed for efficient 3D and 2D representation with 1D separability.-A novel importance sampling strategy for data-driven BRDFs of PDV parameterization.-PDV basis measurement and intuitive BRDF editing from mixture of PDV factors.
We evaluate and compare our models to current stateof-the-art using the measured BRDF data from the MERL database [26], and show that the proposed BRDF models lead to very small approximation errors and produce visually plausible results in renderings.To demonstrate the benefit of our approach, we discuss how the separable representations can be used for applications such as BRDF capture, compact and accurate representations, efficient important sampling of measured BRDFs, and artistic editing of measured materials.

Background
The analysis and new BRDF models presented in this paper build upon a large body of previous work in appearance measurement and modeling.The following subsections describe how our work relates to data-driven BRDF models, existing BRDF parameterizations, and previous approaches for BRDF factorization.

Data-Driven BRDF Models
The way in which light interacts with matter is, in computer graphics, described using the rendering equation [17].In the rendering equation, the scattering at each surface point in a scene is modeled by the BRDF [30].Accurate and efficient BRDF models are essential components in the rendering of complex synthetic scenes for photo-realistic rendering.Measured BRDF data are commonly stored using the standard parameterization which expresses the BRDF directly in spherical coordinates, (ω i , ω o ).Reparameterization of BRDFs has been studied extensively for decades in order to find more efficient and compact representations.A good representation should ideally lead to reducing storage requirement (i.e., provides compactness), an intuitive representation for user editing purposes, and enable efficient importance sampling.Previous works have explored factorization methods in order to find compact representations directly in the standard parameterization [27,44].However, other parameterizations such as the Half-Diff parameterization [39] and the Half-Out parameterization [21] have also been proposed for obtaining more compact representations [7,32,38] and efficient importance sampling procedure [21].However, previous works have not proposed both compact and intuitive representations of BRDF data, allowing for intuitive user edits, flexible and practical measurement setups, and efficient importance sampling.
In this work, we focus on reparameterization of isotropic BRDFs.Our work is inspired by the parametric ABC BRDF model [24], and we compare their projected deviation vector (PDV) parameterization to the more commonly used parameterizations used by Rusinkiewicz [39] and Lawrence et al. [21].The main contributions of this paper are that the resulting BRDF parameterization enables us to represent measured BRDF data with as little as only two 1D factors for adequate representation, and we introduce a new representation of the deviation vector.In contrast to previous data-driven BRDF models, the proposed BRDF parameterization also induces a clear separability in both 2D and 3D representations.Moreover, we investigate factorization methods based on both 2D and 3D representations in order to show that the separability of our parameterization leads to very compact BRDF representations, well approximated with a rank-1 tensor approximation with nonnegative factors intuitive for BRDF editing.
In many applications, it is useful to represent isotropic 3D BRDFs with 2D approximations.2D BRDF representations have shown to be useful both for efficient storage, in acquisition systems and other computer vision problems [13,37].Romeiro et al. [37] proposed a simple transformation from 3D to 2D for measured BRDFs stored in the Half-Diff parameterization.In the Half-Diff parameterization, the parameters (θ h , θ d , φ d ) are reduced to (θ h , θ d ) by taking the average over φ d .Stark et al. [43] presented mathematical relations of (ω i , ω o ) on Barycentric coordinates and proposed three 2D parameterizations.However, these parameterizations are based on unintuitive and complex mathematical functions which limit user intuition.Barla et al. [5] proposed a parameterization which is similar to the PDV; however, there is not, to our knowledge, any BRDF model support for their parameterization.Thus, it would limit the use of its parameterization and importance sampling strategy.

BRDF Factorization
In previous work, it has been shown that measured BRDFs can be factored into basis representations to enable compact representations for efficient rendering and editing.Kautz and McCool [18] proposed the normalized decomposition (ND) method for interactive rendering.In their work, they used a modified Half-Diff parameterization to avoid numerical instability.McCool et al. [27] developed homomorphic factorization on the standard parameterization.The homomorphic factorization is applied on logarithmically transformed data to ensure that the factorized results do not contain negative values.Lawrence et al. [21] proposed the double factorization method on the Half-Out parameterization for material editing applications.The method uses nonnegative matrix Factorization (NMF) [22] on 4D BRDF data.They rearranged 4D BRDFs into 2D BRDFs so that the double factorization could be applied.However, accurate results required a high number of terms in the factorization.Bagher et al. [4] proposed factored BRDF models based on the microfacet BRDF model [11] by using specific weighting functions in optimization process.The resulting models are efficient representations; however, the importance sampling is still an issue on their models; thus, their models require complicated sampling method such as multiple importance sampling (MIS) [40,41,46] or precomputed sampling scheme [25].Ben-Artzi et al. [6] proposed factored BRDF models in the form of linear combinations of wavelet bases to enable real-time BRDF editing and rendering.Their method supports both parametric and analytical BRDF models.An interesting venue for future work is to investigate how our models can be incorporated with their system for editing and rendering.
Tensor factorization methods have also been applied to BRDFs.Sun et al. [44] applied Tucker factorization on the standard parameterization for interactive relighting.Schwenk et al. [42] used CANDECOMP/PARAFAC Decomposition (CPD) method [9,15] on the standard parameterization with an additional dimension which is for wavelength information.The method iteratively applies the rank-1 approximation of the CPD method with repetition of residual factorization.They found that in standard spherical coordinates, the CPD method works better on diffuse and moderately glossy materials; while using the Half-Diff parameterization, CPD works well with glossy materials.A problem, however, is that the method needs at least four factor packs, i.e., the number of iterations on residual factorization, to accurately approximate BRDFs.Later, Ruiters and Klein [38] presented a way to improve optimization of tensor factorizations.They used the Half-Diff parameterization and applied the CPD method on BRDFs transformed to the log-domain using the square of the relative error as the error metric in the optimization.This error metric can provide better approximation compared to using the more common L 2 error metric.However, the output needs up to eight components, i.e., rank-8 approximation, to get accurate results.Bilgili et al. [7] employed the Tucker factorization on 4D BRDFs to represent both isotropic and anisotropic materials in the log-domain.In their work, they simplified the three color channels into one luminance channel and factorized the luminance data.To get RGB BRDF data, a linear regression model is fitted to each of color channels.It is reported that measured materials [26] need up to 13-15 iterations of rank-1 approximations to get accurate results.In contrast to previous work, we show our proposed models can accurately describe BRDFs as a single term rank-1 nonnegative approximation, which allows for both compact and intuitive BRDF models.

BRDF parameterization
The first step in creating a separable BRDF representation is to select a suitable parameterization.In this section, we give an overview of the BRDF parameterizations introduced by Rusinkiewicz [39], Lawrence et al. [7,21], and Löw et al. [24] and analyze their behavior with respect to separability and factorization.We show that they are all good choices for going from 3D to 2D and that the characteristics of the projected deviation vector parameterization [24] make it a particularly good choice for separable BRDF models.This analysis then forms the basis for the new models proposed in this paper.
The values of a BRDF typically exhibit a very high dynamic range [31].To avoid computational problems in such as factorizations and basis representations, it is therefore common to transform the BRDF values to the log-domain.To represent BRDFs efficiently for factorization, we transform the data values using ρ t = log(ρ + 1).

The Half-Difference vector parameterization (Half-Diff)
was introduced by Rusinkiewicz [39], who showed that the specular and retroreflective lobes are aligned with the proposed Half-Diff parameterization.This results in a reduced number of basis coefficients required to accurately represent measured BRDF data and lower memory storage of isotropic BRDFs compared to anisotropic BRDFs.This parameterization has seen widespread use for storing isotropic BRDFs from measured data [26].Figure 1a illustrates the formation of the vectors and the notations used for the Half-Diff parameterization.The parameterization is formed by two vectors, defined in spherical coordinates, the half vector(θ h , φ h ) and the difference vector(θ d , φ d ).The half vector, ω h , is the normalized vector sum of ω i and ω o , given in Eq. 1.The difference vector, ω d , represents the ω i vector in the transformed space.The transformed space is obtained by rotating the hemispherical space so that the half vector is the normal of the transformed space, see Eq. 1.
where rot x,a Y denotes the rotation of Y around x by a • .Moreover, n and b are the normal and binormal vectors to the plane, respectively. Input: Algorithm 1: Standard-to-PDV parameterization conversion.
The Half-Diff parameterization has four parameters which are (θ h , φ h , θ d , φ d ).To represent isotropic BRDFs, φ h can be ignored because isotropic BRDFs are independent of φ h due to symmetry.
The Half-Outgoing vector parameterization (Half-Out) has been used for BRDF factorization in several previous studies [7,21].The Half-Out parameterization is formed by (ω h , ω o ), the half vector, and the outgoing vector.Figure 1b illustrates the vectors and their notations.When ω o is fixed, changing ω h results in moving ω i .Here, the half vector is similar to the Half-Diff parameterization defined by the normalized sum of ω i and ω o , given in Eq. 1.
For isotropic materials, the 4D BRDFs can thus be reduced to 3D BRDFs representations where the BRDFs can be parameterized by either In Sect.3.1, we show illustrations of the Half-Out coordinates and its characteristics.

The projected deviation vector parameterization (PDV)
was used in the parametric ABC BRDF model proposed by Löw et al. [24].The generalized version of the PDV parameterization can be defined by (ω r , D p ) where ω r is represented by spherical parameters (θ r , φ r ) and (d p , φ p ) denote the D p vector.Additionally, ω r is the perfect reflection vector of ω o .The D p vector is formed by the deviation vector between the projections of ω r and ω i on the plane.Figure 1c illustrates the parameter formations and notations used in the PDV parameterization.For isotropic BRDFs, the PDV parameterization has three parameters (θ r , d p , φ p ), where θ r is the zenith angle of the perfect reflection vector from the normal and d p is the length of the D p vector.Its values are [0, 2).Lastly, φ p represents the azimuthal angle of the D p vector in the projection plane.Algorithms 1 and 2 show how the conversions from standard coordinates to the PDV parameterization and vice versa are carried out.
Since d p describes the shape of the BDRF lobe, the quantization over its [0, 2) domain will affect the accuracy.To accurately represent a sharp lobe, very high resolution is required in the quantization.Figure 2 shows plots of two BRDFs demonstrating that most of the variation in the lobes can be found around small d p parameter values, and that  Based on these observations, we use an adaptive quantization.
To find a good quantization that works in the general case, we compute the mean BRDF over a large collection of real materials and compute the nonlinear adapted step length in the quantization using the inversion method, for an overview, see the textbook by Pharr and Humphreys [34].The resulting quantization with 90 nonlinear quantization steps, used in our experiments, is shown in Fig. 3.For the experiments in this paper, we used the MERL BRDF database [26].Although we used all BRDFs in the MERL database, it should be noted that it may be beneficial to group them into classes and use different quantizations for different types of materials.Such analysis, however, is beyond the scope of this paper and left for future work.

Analysis of parameterizations
The different parameterizations described above can all accurately describe isotropic BRDFs using their inherent 3D representation.Since our goal is to reduce the dimensionality of the models required to accurately represent BRDF data, we analyze how well the 3D representation can be represented in 2D.To do this, we first, in Fig. 4, plot the projection of the 3D coordinates on the hemisphere onto the unit disk for regular points in the parameter spaces.The top row shows the Half-Diff parameterization, with parameters (θ h , θ d , φ d ), for a fixed θ h = 30  [24] observed that these circles closely model the iso-contours on the BRDF lobe.This means that the BRDF values along each circle are close to constant.Further conclusions can be drawn from visualizing real BRDF data.Figure 5 illustrates an example from our investigation where the alum-bronze BRDF is plotted in each of the parameterizations, respectively.The three leftmost columns demonstrate how the 2D slices in the 3D representation vary along the third dimension of each parameterization.approximated using nonnegative rank-1 approximations so that the data can be represented as nonnegative univariate functions which lead to intuitive representations for applications such as BRDF capture and editing.
Data factorization can be carried out in several ways.In previous work, nonnegative matrix factorization (NMF) [22] has been a popular choice for factorization of BRDF data, see for example the paper by Lawrence et al. [21].Another commonly used technique is Tucker decomposition [45], which decomposes a tensor into several factors.For 3D data, the approximated data can be computed by a sum of factor products and a core tensor.For the models presented in this paper, we employ a tensor factorization technique called CAN-DECOMP/PARAFAC Decomposition, Canonical Polyadic Decomposition (CPD) [9,15] with the nonnegativity constraint.Compared to the Tucker decomposition, the CPD method is less expensive to compute for random access.Moreover, the Tucker decomposition has a core tensor, [10], which grows exponentially in size if the approximation rank is high.A large tensor core thus leads to higher storage requirements for factorized BRDFs compared to CPD.
For rank-1 approximations, the three techniques are comparable in terms of approximation error.Figure 6 shows a comparison of Tucker decomposition and CPD applied to BRDF data from the MERL database [26] represented using the PDV, Half-Out, and Half-Diff parameterizations for 3D data described in the previous section.Figure 7 shows similar plots but compares NMF and CPD for factorization of 2D data as described in the previous section.Table 1 and 2 show the errors from Figures 6 and 7 in numerical form showing that the PDV parameterization exhibits the smallest mean error and low variance as compared to the Half-Diff and Half-Out parameterizations.Although the error varies between the different factorization models, the approximation accuracy is more or less the same for all three factorization techniques.Given advantages such as a lower computational complexity and efficient random access, we use the CPD method for both 3D and 2D data.Each column shows the BRDF in which one of the parameters was varied.The last column shows the 2D BRDF representations of each parameterization using a higher number of factors for the parameter dimensions, for an overview see [7,19,42].In our models, we have the option of using an iterative method for residual factorization.Since the intuitiveness inherent to the single factor representation is lost, we do not enforce the nonnegative constraint for the residual factorization as this has the tendency to lead to a smaller error.Starting from the original tensor, we constraint the first factorization with the nonnegativity in order to preserve the physical meaning within the BRDFs.For the residual parts, we iteratively factorize the residuals without any constraint, allowing negative values.

New BRDF models
This section presents two new separable BRDF models, in 3D and in 2D, which are based on the above analysis of the PDV parameterization and factorization techniques.The separable models are designed to enable factorization into three or two 1D factors, respectively.

3D separable BRDF model
The first model describes BRDFs using the 3D PDV representation and can be expressed as: where ρ t = log(ρ + 1).F 1,l , F 2,l , F 3,l are the 1D factored functions from the rank-1 iterative CPD method applied to the 3D BRDF representation, and L is the number of iterative terms.As mentioned above, the BRDF can be accurately represented using single factors or a higher number of factors without the nonnegativity constraint to further reduce the approximation error.by using the rank-1 CPD method.This can be formulated as:

2D separable BRDF model
where G 1 , G 2 are the 1D factored functions from the rank-1 CPD method on 2D BRDF representation.

Importance sampling
Efficient BRDF sampling is a key building block in modern Monte Carlo algorithms.For a given BRDF, importance sampling relies on deriving a distribution, or probability density function (PDF) from which samples are drawn during rendering.If the PDF is proportional to the underlying BRDF, the variance of the Monte Carlo estimator is reduced.In this section, we describe how to compute PDF from the PDV factorization.
The goal of computing the PDF is to find the conditional probability, p(d p , φ p |θ o ), so that, for a given viewing direc-     The PDV probability density used for importance sampling can thus be expressed as: To choose the directions, we compute the cumulative distribution function (CDF) as: ( During rendering, we would like to draw sample directions based on a given probability.One of the simplest ways to achieve this is to generate a uniform random variable in range [0,1] and apply the inversion method on the computed CDF.In this case, we generate two uniform random variables (ξ 1 , ξ 2 ) to sample from p φ (φ p ) and p(d p |θ r ).The φ p direction can be computed as: where ξ 1 is the uniform random variable.The magnitude d p = |D p | of the D p vector can be found from the CDF using the inversion method as: where ξ 2 is the second uniform random variable.
It is unavoidable to sometimes choose (d p , φ p ) pairs that result in invalid ω i directions, i.e., falls outside the domain of the parameters.However, when such a (d p , φ p ) pair is chosen, it is straightforward to reject such samples without introducing inaccuracy using rejection sampling, for an overview see the book by Pharr and Humphreys [34].
For clarity, the importance sampling above was described in PDV space.The PDF must thus be transformed to conform  8 A plot of log(J T ) plot is computed by using Eq. 9 and shows the valley caused by the sin(θ i ) factor in the numerator to the Monte Carlo estimator in the standard parameterization using the Jacobian of the mapping: where |J T | is the absolute value of the Jacobian of the transformation.We have p φ (φ p ) = 1 2π because it is uniform random over the φ p parameter.The Jacobian of the mapping is expressed as: where (θ i , φ i ) is the light direction and (θ o , φ o ) is the viewing direction.
Unfortunately, the sin(θ i ) factor in the numerator introduces a valley in the Jacobian surface, see Fig. 8, along which it drops to zeros and exhibits numerical instabilities leading to noise in the renderings.Hence, we propose a smooth version of J T , where the sin(θ i ) factor is omitted: We have found that removing sin(θ i ) does not change the efficiency of the importance sampling significantly and that this approximation leads to high quality rendering results.The effect of the importance sampling is evaluated in the next section.

Results and evaluation
This section presents an evaluation of the new separable BRDF models described in Sect.5, the different parameterizations described in Sect. 3 with respect to their ability to go from 3D to 2D representations, and the effect of the factorization techniques described in Sect. 4 in computing factored models using the different parameterizations.We also show examples of how the models presented in this paper can be used in applications such as storage, efficient BRDF capture and editing.
To measure the error introduced by the models, we use the relative root mean square error (RMS) over all the 100 BRDFs in the MERL database [26].For each combination of BRDF model, parameterization and factorization, we compute the error over the entire hemisphere.In order to generate the same set of sample points and avoid unnecessary interpolation, we first use a random cosine sampling scheme in standard spherical coordinates [34] and then transform each sample into the parameter space of the parameterization where the error is measured.Both the PDV and Half-Diff models are represented at a resolution of 90 × 90 × 360 elements for the (θ h , θ d , φ d ) and (θ r , d p , φ p ) parameter dimensions, respectively.We parameterized the Half-Out parameterization using (θ o , θ h , φ h ) and found that the θ h resolution needs to be higher than 90 for high fidelity renderings.Therefore, for the Half-Out, we used a resolution of 45 × 200 × 360 elements in order to ensure a fair comparison.The supplementary material includes both interactive examples for visual inspection of the proposed models and rendered examples of all BRDFs in the MERL database.
The CPD3D and CPD2D methods were implemented using the N-way Toolbox for MATLAB [2] with a nonnegativity constraint for the first iteration, L = 1.For the experiments, we used the MERL database (including all measured angles).For CPD3D, we used the PDV parameterization at a resolution of 90 × 90 × 180.All the rendering results were produced by using PBRT [34].The fitting time comparison in Table 3 was carried out using an Intel(R) Xeon(R) W-2123 3.60GHz computer with 32 GB memory.Single factor models in 3D and 2D To evaluate the accuracy introduced using the different parameterizations for the single iteration factorization, we examine the separable 3D and 2D models described in Sect. 5 using the PDV, Half-Out, and Half-Diff parameterizations.Figure 6 shows that, in the separable 3D model described in Eq. 2, the PDV parameterization (blue plot) on average leads to a smaller error compared to both Half-Out and Half-Diff for both the Tucker and CPD factorization schemes.The RMS is computed as the error between each 3D representation and its corresponding factored model.Visual inspection of the rendered BRDFs reveals that the PDV parameterization produces significantly better results compared to the other parameterizations, but the Half-Diff parameterization has a tendency to produce better results for diffuse materials.An example of this is shown in the second columns (CPD3D-Iteration-1) in Fig. 10, where renderings of the chrome and special-walnut-224 BRDFs are visualized for different combinations of parameterizations, model dimensionality and number of iterations.Figure 7 shows the errors introduced when going from the 3D representations of the PDV, Half-Out, and Half-Diff parameterizations to the corresponding factored models in 2D described in Eq. 3. Similarly to the factored 3D models, it is evident that the PDV parameterization on average leads to a smaller error compared to Half-Out, and Half-Diff.Visual inspection, see the sixth columns (CPD2D) in Fig. 10, shows that the PDV parameterization, also in the 2D case, tends to handle glossy materials with small errors and that the Half-Diff parameterization is better at diffuse materials.In summary, the evaluations and the renderings show that the proposed factored 3D and 2D models produce both small numerical errors and visually plausible rendering results using only a single factor rank-1 approximation, and that PDV and Half-Diff are suitable parameterizations for glossy and diffuse materials, respectively.factors for higher accuracy.As described in Sect.4, this can be carried out using iterative factorization of the residual error so that each factor is represented using L basis vectors.

Iterative factor 3D model
Comparing the error plots in Fig. 6b using L = 1 iteration to that in Fig. 9 using L = 10 iterations shows that both the maximum and mean errors are reduced significantly for all parameterizations, please see Table 4.It is also evident that the different parameterization on average starts to perform on par with each other in terms of numerical error.This is likely due to that the increased expressive power of the higher-order models can compensate better for the residual errors.The fourth column (CPD3D-Iteration-10) in Fig. 10 reveals similar behavior for the different parameterizations as for the L = 1 iteration renderings in the second column.Finally, Fig. 11 shows examples of how the accuracy varies with an increasing number of L = 10, 15, and 20 iterations.Comparing to the models introduced by Lawrence et al. [21] and Bilgili et al. [7], the color-coded insets and the peak signalto-noise ratio (PSNR) [36] measures show that our iterative model produces state-of-the-art results.The different models lead to different storage requirements.Table 5 compares the storage requirements for each of the models used in Fig. 11 as well as the reference measured BRDF data.Table 6 and Figure 12 compare the proposed 3D models to the Naive model proposed by Bagher et al. [4].The blue curve in Figure 12 shows the error plot of Naive model sorted in ascending order.We chose the Naive model as a representative method from Bagher et al. since it is the simplest method among their models and due to the fact that Bagher et al. models perform similarly to each other.Our 3D model with L = 15 performs on par with the Naive model.When L is increased, our model performs overall better than to the Naive model.Table 6 summarizes the mean, minimum and maximum errors computed over all BRDFs in the MERL database comparing our 3D models to Bagher et al. for different numbers of iterations, L.
Importance sampling Figure 13 compares the importance sampling technique developed for our factored models described in Sect.5.1 to the techniques presented by Lawrence et al. [21], Edwards et al. [12], and Bilgili et al. [7].The Figure shows renderings of the Princeton scene generated using the PBRT rendering system, described in the book by Pharr and Humphreys [34], with 256 path samples per pixel.The scene consists of Nickel, Yellow-matte-plastic, and Blue-metallic-paint materials.The rendered images show that our importance sampling technique performs better than Lawrence et al. [21] and Edwards et al. [12] techniques numerically, but our technique gives numerically worse performance than Bilgili et al. [7] technique.A natural next step, however, would be to derive a more accurate approximation of the Jacobian which potentially would lead to even better results.The color shift on and near the teapot in Fig. 13c is due to Bilgili et al.'s representation inaccuracy.Because this representation factorizes intensity channel and then estimates color values linearly as a sub-procedure by a few linear coefficients.Table 7 shows the rendering time for the images in Fig. 11 (in seconds).Even with our higher accurate models (L > 15), our factored models have comparable computational cost compared to Edwards et al. [12], better computational cost compared to the technique presented by Bilgili et al. [7], but a higher computational cost compared to the technique presented by Lawrence et al. [21].Note that our   factored models compute probabilities on the fly and do not rely on precomputed lookup tables for importance sampling, but many previous BRDF representations (i.e., Bilgili et al. [7] and Lawrence et al. [21]) rely on precomputed lookup tables for importance sampling.Figure 14 shows the importance sampling efficiency evaluated using the logarithmic MSE (mean squared error) from a sphere rendered in a diffuse (constant) white environment.The results show importance sampling efficiency of various techniques, which are similar to the results demonstrated in Fig. 13.
Limitations Our approach is dependent on the quantization of the D p vector.In the experiments presented here, we computed the d p step length based on the average of all materials in the MERL database.We have noticed that this results in unsuitable quantization for some materials, see special-walnut-224 in Fig. 10.This could be alleviated by computing different nonlinear quantizations for different material classes.Additionally, in order to make a comparison with Bagher et al. [4], we included all angles of the materials in the MERL database.Previous studies (e.g., Ngan et al. [28] and Löw et al. [24]) have pointed out that the measured values in the MERL database are increasingly unreliable toward grazing angles.It would therefore be interesting to evaluate the accuracy of the models within a restricted range of angles.Moreover, the Jacobian term used in the importance sampling is an approximation to the exact solution due to numerical instabilities.This could be improved by deriving a better approximation.

Applications
This section shows examples of how the PDV factored models can be used for efficient measurement of isotropic BRDFs and editing of measured BRDF through the individual PDV factors.

BRDF reconstruction from factor measurements
Using the factored 2D model described in Eq. 3, it is possible to capture a full isotropic BRDF by measuring the scaling of the lobe described by G 1 (θ r ) and the shape of the lobe described by G 2 (d p ) in a planar slice of the BRDF.In the illustration of the 2D PDV parameterization in Fig. 5l, this would correspond to measure the 1D G 2 (d p ) factor representing the vertical direction and the 1D G 1 (θ r ) factor representing the horizontal direction in the 2D matrix.As illustrated in Fig. 15 (left), the horizontal factor G 1 (θ r ) can be measured by moving a sensor and a calibrated light source in opposite  7 The rendering time (in seconds) for the renderings shown in Fig. 11 and a  The Princeton scenes were rendered using 256 samples per pixel.Each single-material scene comprises a sphere rendered using 1024 samples per pixel at a resolution of 256 by 256 pixels.All renderings and timings were carried out using the same computer environment as described in Sect.6  factor needs to be normalized before it can be modulated with G 1 (θ r ) to form the full 2D matrix describing the BRDF.
Figure 16 shows the reconstruction errors (red plot) obtained by simulating capture and reconstruction as described above for all BRDFs in the MERL database, and as a reference the

Editing of measured BRDFs
The single iteration (L = 1) nonnegative factor models in 2D and 3D enable intuitive editing.For example, in the 2D model described in Equation 3, G 2 (d p ) can be edited to change the shape of the lobe, and the G 1 (θ r ) can be edited to change effects such as graz-

Conclusions and future work
This paper presented a study of three different BRDF parameterizations and proposed two new factored BRDF models for 3D and 2D representations of isotropic BRDFs.The study showed that the Half-Diff [39], Half-Out [21], and PDV [24] parameterizations are all suitable for 2D representation of isotropic BRDFs with only small approximation errors.It also showed that the PDV parameterization structures the data in a way suitable for factorization.The evaluation demonstrated that the new models are flexible in terms of which parameterization is used and that they can represent measured BRDF with low numerical approximation errors, and produce visually plausible renderings.For efficient rendering, the paper also presented a new Monte Carlo importance sampling scheme based on the factored models.The evaluation and the results point toward a number of interesting venues for future work.One interesting area is to further study how the choice of parameter quantization affects the accuracy of the factorization.Another important direction is to investigate how low-dimensional factored models can be adapted for accurate representation of anisotropic BRDFs.Finally, we will use the factored models to develop efficient BRDF capture systems and investigate how they can be incorporated in computer vision applications where BRDFs are needed to be characterized in the wild.

Fig. 2 Fig. 3
Fig. 2 The figures contain examples of BRDF plots in the PDV parameterization.Both BRDFs are from a fixed angle of θ r = 45 • and φ p = 0 • .Vertical axis is BRDFs scaled by logarithmic function and Horizontal axis is d p resolution.The left figure shows the BRDF of alum-bronze which represents the class of glossy materials.The right figure shows the BRDF of blue rubber which represents the class of diffuse materials.It is apparent that BRDFs in the PDV parameterization aligned mostly around small d p Figure 5a-d illustrates the BRDF data in the Half-Diff parameterization along the φ d , Fig. 5e-h illustrates the BRDF data in the Half-Out parameterization along the φ h dimension, and Fig. 5i-l illustrates the BRDF data in the PDV parameterization along the φ p dimension.The right columnshows the corresponding 2D representation computed as the mean over all slices in the third parameter dimension for each parameterization.The plots illustrate that for all three parameterizations, the BRDF values are stable along the third dimension, which supports the conclusion that they are suitable choices for going from a 3D to a 2D representation.Visual inspection of the BRDFs in the MERL database also shows that data represented in the PDV parameterization is better aligned (horizontally and vertically) in the 2D plane compared to the Half-Diff and Half-Out parameterizations.As we show in the next section, this structure makes it possible to accurately factorize the 2D representation into a 1D + 1D factored BRDF model.

Fig. 4
Fig. 4 Figure a-d illustrates the coordinate samples of the Half-Diff parameterization with fixed angles of θ h = 30 • , 70 • and varying θ d , 2 degrees apart, and rotating φ d ∈ (0, 2π).Figure a and b were sampled on the hemisphere.Figure c and d were sampled on the disk.Figure e-h illustrates the coordinate samples of the Half-Out parameterization Fig. 4 Figure a-d illustrates the coordinate samples of the Half-Diff parameterization with fixed angles of θ h = 30 • , 70 • and varying θ d , 2 degrees apart, and rotating φ d ∈ (0, 2π).Figure a and b were sampled on the hemisphere.Figure c and d were sampled on the disk.Figure e-h illustrates the coordinate samples of the Half-Out parameterization

Fig. 6
Fig.6 Both plots show errors between the BRDF data and its rank-1 approximation.The Tucker decomposition was used in (a).The CPD method was used in (b).Both methods were applied on three different parameterizations.Blue line represents errors from the PDV parameterization.While red and yellow lines represent errors from the Half-Out and the Half-Diff parameterizations, respectively

Fig. 7
Fig.7 Both plots show errors between the 3D BRDF data and its factored 2D representations.The NMF was used in (a).The CPD method was used in (b).Both methods were applied on three different parameterizations.Blue line, red line, and yellow line represent errors from the PDV, the Half-Out, and the Half-Diff parameterizations, respectively Fig.8A plot of log(J T ) plot is computed by using Eq. 9 and shows the valley caused by the sin(θ i ) factor in the numerator

Fig. 11 A
Fig. 11 A comparison in which the Princeton scene is rendered using b, c, d and e our 3D model with the PDV parameterization with L = 5, L = 10, L = 15, and L = 20 iterations, respectively, and f, g and h the Edwards et al., the Lawrence et al. and the Bilgili et al.BRDF models.All images were rendered at 262144 samples/pixel using a path tracing algorithm.Insets show a scaled color-coded difference between the reference image and the rendered image for better comparison.For higher disparity, the color-coded difference images are scaled by 5. Below each image we also report the PSNR value

Fig. 12
Fig. 12 To compare the quality of our models against Bagher et al. model, we show the error plots of our models with L = 5, 15, 20, 25 and the error of Bagher et al.'s naive model.Here it is sorted in ascending order of the naive model errors, blue line

Fig. 13
Fig. 13 The Princeton scene rendered for visual comparisons of sampling efficiency.a, b, c, and d were rendered using the Edwards et al., the Lawrence et al., the Bilgili et al., and our factored BRDF models, respectively.All images were rendered at 256 samples/pixel by using

Fig. 14
Fig. 14 The MSE plots of Nickel BRDF in white constant environment demonstrate the importance sampling efficiency of various methods

Fig. 15
Fig. 15 The factors G 1 (θ r ) and G 2 (d p ) of our separable 2D model can be captured

Fig. 16 (Fig. 17
Fig. 16 (red) The reconstruction error for all BRDFs in the MERL database and (blue) as a reference the approximation error of the factored 3D model using L = 1 iteration Fig. 17Three materials were rendered to compare between the original BRDF data and its reconstruction

Fig. 18
Fig.18 The single iteration 2D model 1D factors G 1 and G 2 from the left and middle columns are mixed to create new material appearances • , 70 • .The plots show how ω o moves on the hemisphere and in the projection plane when ω d is varied.The middle row shows the Half-Out parameterization, with parameters (θ o , θ h , φ h ), and a fixed θ o = 30 • , 70 • .The plots show how ω i changes as a function of ω h .Finally, the bottom row shows the same plots for the PDV parameterization, with parameters (θ r , d p , φ p ).By varying d p and φ p with a fixed angle of ω r = 30 • , 70 • , only ω o is changed and forms circles in the projection plane.In fact, for smooth surfaces Löw et al.

Table 1
The table shows mean, variance, minimum and maximum values of log relative errors of all 100 materials in the MERL database for each method, showing in Fig.6, θ o and θ r can be interchanged.We can thus rewrite the conditional probability as: p(d p , φ p |θ r ).Since the φ p factor is constant up to the truncation of the unit circle, the condi- tional probability can be factored into the multiplication of two PDFs, p(d p |θ r ) p φ (φ p ), where p φ (φ p ) is a uniform PDF.

Table 2
The table shows mean, variance, minimum and maximum values of log relative errors of all 100 materials in the MERL database for each method, showing in Fig.7is the resolution of the d p quantization, in our case 90, and Δd p j is the d p step size.Note that the d p step size is nonlinear due to the quantization described in Sect.3.
The term ρ(θ r , d p , φ p ) can be computed by using Eq. 2 when L = 1.

Table 5
The storage requirements for the BRDF models compared in Fig.11