This document accompanies the SIGGRAPH
2004 paper that describes a factored BRDF representation useful
for efficient importance samplingthe notation used here is the same
as in the paper. This document contains the following information:

Factored BRDF File Format 
This file contains the values of the terms in the factored BRDF representation computed using the NMF algorithm. 
The header lists the parameterization of the incoming hemisphere, the resolution of the original data matrix (and, consequently, the resulting factors) and the number of terms used during the factorization. 
The intensity data gives the actual values of each 'pixel' in each factored term to form an approximation of the intensity of the original BRDF over its entire domain (intensity = (red + green + blue)/3.0). The resolution of each term will be equal to the respective resolution in the original data matrix (i.e. it will be the same as those values listed in the header). 
The color data gives the values of each 'pixel' in the factors required to account for the color of the BRDF. Although not described in the paper, the details of how wavelength dependence of the BRDF is handled is covered in this document. 
HEADER  
param: [type of parameterization]  In most cases, this is either 'halfangle' or 'spherical' 
outer terms ('J' in paper): [# of outer terms]  The number of terms in the first factorizationsee paper, Section 4 for details 
inner terms ('K' in paper): [# of inner terms]  The number of terms in the second factorizationsee paper, Section 4 for details 
The total number of terms, 'L' in the paper will be #outer_terms X #inner_terms  
theta_out_count: [# samples along theta_out]  The number of samples along the elevation angle of the outgoing direction 
phi_out_count: [# samples along phi_out]  The number of samples along the azimuthal angle of the outgoing direction 
theta_p_count: [# samples along theta_p]  The number of samples along the elevation angle of the incomging direction for the current parameterization 
For example, if param = halfangle, this is the number of samples along theta_h  
... , if param = spherical, this is the number of samples along theta_i  
z_count: [# samples along z]  The number of samples in zdimension, where z = 1cos(theta_p) for the given parameterizationthis only shows up in sampling 
phi_p_count: [# samples along phi_p]  The number of samples along the azimuthal angle of the incoming direction for the current parameterization 
INTENSITY DATA  
F_0(theta_o,phi_o)  The 2D factor dependent on the view only, Equation 1 in the paper (stored so that columns/rows vary with phi_out/theta_out respectively). 
G_0(theta_p,phi_p)  The 2D factor dependent on the incoming direction for the given parameterization (stored so that columns/rows vary with phi_p/theta_p respectively). 
u_0_0(theta_p)  The 1D factor dependent on the elevation angle of the incoming direction, see Eq. 3 in the paper 
v_0_0(phi_p)  The 1D factor dependent on the azimuthal angle of the incoming direction, see Eq. 3 in the paperNOTE: only intensity values) 
.  
.  
.  
u_0_{K1}  
v_0_{K1}  
F_1(theta_o,phi_o)  
G_1(theta_o,phi_o)  
u_1_0  
v_1_0  
.  
.  
.  
u_1_{K1}  
v_1_{K1}  
.  
.  
.  
F_{J1}(theta_o,phi_o)  
G_{J1}(theta_o,phi_o)  
u_{J1}_0  
v_{J1}_0  
.  
.  
.  
u_{J1}_{K1}  
v_{J1}_{K1}  
COLOR DATA  
Fc(theta_o,phi_o)  The 2D color factor dependent on the view only, see below for details (stored so that columns/rows vary with phi_out/theta_out respectively and each entry contains RGB written in that order). 
Gc(theta_p,phi_p)  The 2D color factor dependent on the incoming direction only, see below for details (stored so that columns/rows vary with phi_out/theta_out respectively and each entry contains RGB written in that order). 
EOF 
CDF File Format 
This file contains the values of the numerical Cumulative Distribution Functions (CDFs) derived from each term in a factored BRDF representation. These terms can be used to generate samples of the incoming hemisphere proportional to the shape of the factored BRDF for efficient variance reduction in standard global illumination algorithms. 
CDF DATA  
P(phi_pterm 1 selected)  This is the CDF over the incoming azimuthal angle if term 1 is selected 
P(zterm 1 selected)  This is the CDF over the latitude, z, if term 1 is selected 
.  
.  
.  
P(phi_pterm L selected)  Recall there are a total of L=J*K terms in the entire factorization 
P(zterm L selected) 
BRDF Reconstruction 
Assume we want to reconstruct the value of the BRDF at some incoming/outgoing direction (i.e. omega_in/omega_out). We first compute their elevation and azimuthal values for the given parameterization of the hemisphere. The outgoing hemisphere is always parameterized by the traditional spherical parameterization (i.e. theta_out and phi_out are w.r.t. to the local surface normal). However, the parameterization over the incoming hemisphere can vary. To demonstrate, assume the incoming direction is parameterized w.r.t. the halfangle vector. In this case: omega_h = (omega_in + omega_out) / norm(omega_in + omega_out). We then compute the elevation and azimuthal angles, theta_h and phi_h, for this omega_h vector. Now we have the necessary values of theta_out, phi_out, theta_p (=theta_half) and phi_p (=phi_half) to reconstruct the value of the BRDF. 
Next we must transform these angles into indices into each factored term. This will depend on the resolution at which each dimension is sampled (i.e. the respective values of theta_out_count, phi_out_count, theta_p_count, phi_p_count listed above). Now we are ready to compute the value of each tabulated factor at this incoming/outgoing pair to reconstruct the intensity of the BRDF: 
BRDF_intensity(omega_in,omega_out) = SUM(j=1toJ)Fj(theta_out,phi_out)Gj(theta_p,phi_p) 
This actually evaluates the intensity of the BRDF*cos(theta_out). We factor this original form for sampling reasons. 
To account for color in the BRDF we observe the following
relationships:

To this end, we compute a singleterm factorization of the quantities BRDF_red/BRDF_intensity, BRDF_grn/BRDF_intensity, BRDF_blu/BRDF_intensity using the NMF algorithm at the same resolution of the original BRDF (NOTE: we use our current approximation for the BRDF's intensity). These 2D factored terms are denoted Fc(omega_out) and Gc(omega_in) above. To compute the red value of the BRDF at this omega_in/omega_out pair, you do the following: 
BRDF_red(omega_in,omega_out) = Fc(theta_out,phi_out).Red()*Gc(theta_p,phi_p).Red()*BRDF_intensity(omega_in,omega_out) 
The value of the other 2 color channels are computed analogously. 
Sampling 
We are interested in generating incoming directions for a known view (i.e. some value of omega_out). To do this, we first evaluate each viewdependent factor at this omega_out (this is done for each "outer" term in the factorization): 
f[1] = F_1(theta_out,phi_out) 
f[2] = F_2(theta_out,phi_out) 
. 
. 
. 
f[L] = F_L(theta_out,phi_out) 
We then select a term randomly with probability equal to the relative 'energy' in each term. For example, the i^th term is selected with probability equal t= f[i]/(SUM(j=1toL)f[j]). 
Once the term is selected we use the corresponding 1D CDFs in
both phi_p and z (z = 1cos(theta_p)) to generate an incoming
direction (i.e. we use P(phi_pi) and P(zi) as described above. Specifically, we generate two
uniformly distributed random variables (zeta_1 and zeta_2) and use a
binary search to invert the value of each numerical CDF at these
values:

We invert the measure of latitude to elevation (this change of variable accounts for the fact that domega = sin(theta)dthetadphi). We obtain theta_p = acos(1z). Using this value of phi_p and theta_p we can generate a direction w.r.t. the current parameterization of incoming directions. For the halfangle parameterization these values correspond to phi_h and theta_h respectively. Using this halfangle direction and the current view we can generate an incoming direction, omega_in. 
To compute the probability of this sample, we must sum the marginal probabilities over all the terms in the factorization. Intuitively speaking, we must sum over each term because this particular sample (omega_in) had a nonzero chance of being generated by any of term in the factorization. Specifically: 
Prob(omega_in) = SUM(i=1toL)[(f[i]/(SUM(j=1toL)f[j])) * Prob(phi_pi) * Prob(zi) 
The probabilities, Prob(phi_pi) and Prob(zi), are equal to the slope of the line segment of the numerical CDF at the corresponding position of phi_p and z in their domains and can be computed in constant time. 
Lastly, we must account for the effect of reparameterizing the incoming hemisphere (i.e. the measure domega_p does not always equal domega_in over the incoming hemisphere when the mapping between the two is nonlinar). For the case of the halfangle parameterization, we scale the probability by 1/(4*dot(omega_in,omega_half)). See the paper for details. 