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{
public:
float x, y, z;
float theta, phi;
};
MontSample montSamples[montSampleSqrtNum * montSampleSqrtNum];
//local coordinate samples
void initMontSamples(void){
srand((int)time(0));
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 +
xi[2]*montSamples[i].z));
}
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;
else
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.
My Research Projects 