How I built the diffuse component lookup table for Kajiya-Kay model (Updated!)

The environment lighting was approximated by a series of SRBF lights. I built the diffuse component lookup table, and then rendered the model shaded by one SRBF light through texture fetching in the shader. And compared with the result of direct computation using Kajiya-kay model to verify the lookup table.

Following the description in the paper Interactive Hair Rendering Under Environment Lighting, the bidirectional scattering function for Kajiya-Kay model is

S(ωi, ωo)=Kd + Kscosp(θi+θo)/cosθi

The diffuse component of the convolution integral

is independent of ωo and the azimuth angle φ of the SRBF center ξ, and thus can be precomputed as a 2D table Id(cosθξ, 1/λ).

Because this integral integrates over the unit sphere, in order to compute it, I first generated a lot of sample locations on the unit sphere. But before I had done that, I redefined the coordinate system for the coordinate system used in the original paper is a little misleading.

The coordinate system used in the original paper is

Personally, it is better to assume that the u axis corresponds to the z axis in the regular spherical coordinate system we all familiar with, and w corresponds to y, v corresponds to x.

Here is the regular spherical coordinate system. φ ranges from 0 to 2π. θ ranges from -π/2 to π/2 and is positive in the upper hemisphere (θ is the angle between OP and its projection onto the X-Y plane).

The following code were used to generate sample locations over the unit sphere

const int montSampleSqrtNum = 50;//the total number of samples is 50*50

class MontSample{


float x, y, z;

float theta, phi;


MontSample montSamples[montSampleSqrtNum * montSampleSqrtNum];

//local coordinate samples

void initMontSamples(void){


for(int i=0; i<montSampleSqrtNum; ++i)

for(int j=0; j<montSampleSqrtNum; ++j)


float x = (i+rnd(1.0))/montSampleSqrtNum;

float y = (j+rnd(1.0))/montSampleSqrtNum;

float theta = 2.0 * acos(sqrt(1.0 – x)) – PI/2;//acos(sqrt(1.0-x)) returned a value in the region [0, π/2], after multiplying by 2.0 and then subtracting π/2, the result lay in [-π/2, π/2]

float phi = 2.0 * PI * y;

montSamples[j+i*montSampleSqrtNum].theta = theta;

montSamples[j+i*montSampleSqrtNum].phi = phi;

montSamples[j+i*montSampleSqrtNum].x = cos(theta) * cos(phi);

montSamples[j+i*montSampleSqrtNum].y = cos(theta) * sin(phi);

montSamples[j+i*montSampleSqrtNum].z = sin(theta);



Here comes the code to build the lookup table.

void i_dBuild(void){

float phi = 0;//the diffuse component is independent of the azimuth angle φ of the SRBF center, so any value can be used here.

for(int h=0; h<lambdaRes; ++h)//resolution for the bandwidth λ is 64

for(int w=0; w<cosRes; ++w)//resolution for cosθξ  is 32.


float integral = 0;

float cosTheXi = -1.0f + (w*2.0f)/cosRes;// w÷cosRes only results in [0, 1]. Here it has been adjusted to fit into [-1, 1], the range for cosine.

float theXi = acos(cosTheXi) – PI/2;//acos() returns [0, π], subtract π/2 to fit into [-π/2, π/2].

float lambda = 64.0 * lambdaRes / (3.0*h + lambdaRes);

for(int i=0; i<montSampleSqrtNum*montSampleSqrtNum; ++i)


float xi[3] = {cos(theXi)*cos(phi), cos(theXi)*sin(phi), sin(theXi)};

integral += cos(montSamples[i].theta) * exp(-lambda)*exp(lambda*(xi[0]*montSamples[i].x + xi[1]*montSamples[i].y +



i_d[w+h*cosRes] = integral;



The original paper used 1/λ  as a lookup argument. Because λ lies in [16, 64], so 1/λ  lies in [1/64, 1/16], thus the following equation can map 1/λ  to [0, 1], which is the required range for a lookup argument:

(1/λ -1/64)*64/3

Therefore, in the above code, λ can be solved for 64.0 * lambdaRes / (3.0*h + lambdaRes);

In the fragment shader, the code used to fetch the texture value is

vec2 uvId;

uvId.y = (1.0/lambda – 1/64.0)*64.0/3.0;

if(dot(l, t)>0)//t is the normalized tangent, l is the light direction.
    uvId.x = sin(acos(dot(l, t)))/2.0 + 0.5;
    uvId.x = -sin(acos(dot(l, t)))/2.0 + 0.5;

float diffuse = texture(i_dTex, uvId).r;

Run the program with only the diffuse component activated, I got the following result

The reference result I got from directly computing using Kajiya-Kay model

The first image shows false location of the low intensity area.

And the code for the diffuse component of Kajiya-Kay model:

diffuse = diffCoeff * sqrt(1 – dot(t, l) * dot(t, l));

BTW: I did not see the gain of performance by using precomputed tables in this case.


The original paper defines θ to be in the interval [-π/2, π/2], which is monotonic for sine but not cosine, therefore, if use cosθ to index into the lookup table as well as generate the lookup table ahead, the information of whether the angle is positive or negative will be lost.

So simply replace with the following two lines in the corresponding place in i_dBuild()

float sinTheXi = -1.0f + (w*2.0f)/cosRes;
float theXi = asin(sinTheXi);

and replace with

uvId.x = dot(l, t)/2.0 + 0.5;

in kajiyaShader.fp.

The rendering result will be correct this time.


Analysis of Marschner’s model

Here is my understand regarding the azimuthal component of Marschner’s model. All the analysis are carried on a cross section of the hair.

I posted two figures for later reference. Figure 1 shows the scattering geometry, and Figure 2 the cross section.

Figure 1: Scattering geometry

Figure 2: Geometry for scattering from a circular cross section

The goal of Marschner’s model is to solve for the outgoing intensity. Please refer to Figure 2 first, the incoming ray has been scattered into three ways, some part of the energy has been reflected directly; some enters the hair, and leaves at the other side of the surface; some enters the hair, reflected at the inner surface and final refracted out. Likewise, when given the outgoing angle phi(the relative azimuth phi_r – phi_i), we can reverse this process to find out the directions of light that contribute to the outgoing intensity. For instance, in Figure 3, the intensity of the outgoing direction (in blue) is the contribution from three (not necessary to be exactly three)incoming beams (in  red).

Figure 3: Light from three paths that contribute to one outgoing direction

Back to the second figure, setting aside attenuation for the moment, power from a small interval dh in the incident beam is scattered into an angular interval d/phi in the exitant intensity distribution


To complete this equation, we have to take into account attenuation caused by absorption and reflection by introduce an attenuation factor A(p, h) in front of the intensity contributed by a path.

where p = 0 stands for surface reflection R, 1 for refractive transmission TT, and 2 for internal reflection TRT. Refer to Figure 3 for a visualization of the three paths.

Consider Path 1 in Figure 3, the light reflects directly, so the Fresnel factor can be applied to account for the reflected energy, thus A(0, h) = F(eita, gamma_i);

For Path 2, the light first undergoes refraction, thus the energy enters the hair should be (1-F(eita, gamma_i)) multiple the original energy. Then it travels a distance of 2*cos(gamma_t) (with the radius of the hair being unit length), assume sigma_a to be the volumn absorption per unit length, the absorption factor T should be exp(…), finally, it refracts out of hair, so (1-F(eita, gamma_t)) is used to account for the refracted energy.

Path 3 is similar to Path 2, except that the light reflected once at the inner surface of hair, and travels two times of the distance as it in Path 2. Personally, I would write A(2, h) as

rather than

given by the paper, where p = 2.

One more thing to mention. Before we carry on analysis, we have to project the scattering geometry from 3D to 2D, as well as various parameters, including the index of fraction, and each internal segment should be lengthened by a factor of 1/cos(theta_t). 

 The End