• 9
• 10
• 10
• 11
• 16

# Spherical Harmonics - analytical solution

This topic is 1059 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

So, I've gone through this paper (http://www.inf.ufrgs.br/~oliveira/pubs_files/Slomp_Oliveira_Patricio-Tutorial-PRT.pdf). It's all perfectly clear how everything works (or at least I think so), I implemented the stuff - it works. The only thing that bugged was that I had to do pretty lot of computation to just project some analytical lights (as the code in the paper shows you integrate the light function through a number of samples, which could be a lot) and I was pretty sure I had seen SH to be "simpler". And I found chapter 2.15 in ShaderX3. This is the code that author uses for projecting a single directional light into SH:

static void CalculateCoefficentsPerLight(LightStruct *light)

{

SHCoeff[0].red += light->colour[0] * fConst1;

SHCoeff[0].green += light->colour[1] * fConst1;

SHCoeff[0].blue += light->colour[2] * fConst1;

SHCoeff[1].red += light->colour[0] * fConst2 * light->direction[0];

SHCoeff[1].green += light->colour[1] * fConst2 * light->direction[0];

SHCoeff[1].blue += light->colour[2] * fConst2 * light->direction[0];

SHCoeff[2].red += light->colour[0] * fConst2 * light->direction[1];

SHCoeff[2].green += light->colour[1] * fConst2 * light->direction[1];

SHCoeff[2].blue += light->colour[2] * fConst2 * light->direction[1];

SHCoeff[3].red += light->colour[0] * fConst2 * light->direction[2];

SHCoeff[3].green += light->colour[1] * fConst2 * light->direction[2];

SHCoeff[3].blue += light->colour[2] * fConst2 * light->direction[2];

SHCoeff[4].red += light->colour[0] * fConst3 * (light->direction[0] * light->direction[2]);

SHCoeff[4].green += light->colour[1] * fConst3 * (light->direction[0] * light->direction[2]);

SHCoeff[4].blue += light->colour[2] * fConst3 * (light->direction[0] * light->direction[2]);

SHCoeff[5].red += light->colour[0] * fConst3 * (light->direction[2] * light->direction[1]);

SHCoeff[5].green += light->colour[1] * fConst3 * (light->direction[2] * light->direction[1]);

SHCoeff[5].blue += light->colour[2] * fConst3 * (light->direction[2] * light->direction[1]);

SHCoeff[6].red += light->colour[0] * fConst3 * (light->direction[1] * light->direction[0]);

SHCoeff[6].green += light->colour[1] * fConst3 * (light->direction[1] * light->direction[0]);

SHCoeff[6].blue += light->colour[2] * fConst3 * (light->direction[1] * light->direction[0]);

SHCoeff[7].red += light->colour[0] * fConst4 * (3.0f * light->direction[2] * light->direction[2] - 1.0f);

SHCoeff[7].green += light->colour[1] * fConst4 * (3.0f * light->direction[2] * light->direction[2] - 1.0f);

SHCoeff[7].blue += light->colour[2] * fConst4 * (3.0f * light->direction[2] * light->direction[2] - 1.0f);

SHCoeff[8].red += light->colour[0] * fConst5 * (light->direction[0] * light->direction[0] - light->direction[1] * light->direction[1]);

SHCoeff[8].green += light->colour[1] * fConst5 * (light->direction[0] * light->direction[0] - light->direction[1] * light->direction[1]);

SHCoeff[8].blue += light->colour[2] * fConst5 * (light->direction[0] * light->direction[0] - light->direction[1] * light->direction[1]);

}


Later on, these nine SHCoeffs are passed to the vertex shader to compute the final illumination:

     vec3 norm = gl_Normal;

color  = Coefficients[0];
color += Coefficients[1] * norm.x;
color += Coefficients[2] * norm.y;
color += Coefficients[3] * norm.z;
color += Coefficients[4] * norm.x * norm.z;
color += Coefficients[5] * norm.y * norm.z;
color += Coefficients[6] * norm.x * norm.y;
color += Coefficients[7] * (3.0 * norm.z * norm.z  - 1.0);
color += Coefficients[8] * (norm.x * norm.x  - norm.y * norm.y);


This stuff looks much simpler and I don't see any need for sophisticated integration, nor even any dot(lightVector, normal) term (which is present in the PRT Tutorial paper in the projection of the transfer function). So my question is simply about the derivation of the ShaderX3 article's constants and how this relates to the "general" solution presented in the PRT Tutorial paper?

##### Share on other sites

nor even any dot(lightVector, normal) term

There is definitely a dot product between those two terms there, it's just split between the precalculation and the shader.

Recall that dot(L, N) = L.x*N.x + L.y*N.y + L.z*N.z.

The light direction is stored in slots 1, 2, and 3 of the pre-calculated Coefficients array, which are each multiplied by the matching component of the surface normal in the shader, and all are finally added together...

##### Share on other sites

Oliveira uses numerical integration, because in his tutorial he projects a light probe and shadowed lights into SH basis. Both usually don't have an analytical form, so he can't leverage orthogonality of SH basis and solve that analytically.

##### Share on other sites

The light direction is stored in slots 1, 2, and 3 of the pre-calculated Coefficients array, which are each multiplied by the matching component of the surface normal in the shader, and all are finally added together...

Correct :).

Oliveira uses numerical integration, because in his tutorial he projects a light probe and shadowed lights into SH basis. Both usually don't have an analytical form, so he can't leverage orthogonality of SH basis and solve that analytically.

Right. So the question now is what my light and transfer functions should look like so after analytical integration I end up with coefficients from the ShaderX3 article?

##### Share on other sites

For Lambertian BRDF transfer function is cos(theta) and directional light is projected by evaluating the SH basis function in the light's direction and scaling results by light's color. Those coefficients from the ShaderX3 article are derived from the spherical harmonic basis functions (Yml).

Edited by Krzysztof Narkowicz

##### Share on other sites

For Lambertian BRDF transfer function is cos(theta) and directional light is projected by evaluating the SH basis function in the light's direction and scaling results by light's color. Those coefficients from the ShaderX3 article are derived from the spherical harmonic basis functions (Yml).

Cos(theta) happens to have a circular invariance. With correct basis (ie where z is an axis of rotation invariance) the SH is a zonal spherical harmonics ie all its Yml coefficients are 0.f if m != 0.

The SH representation of a cosfunction for a light whose direction is toward z is thus the sum of all Ym0 for m >= 0, scaled by light's color.

Now not all lights are directed toward z but thanks to rotation invariance property there is a "simple" formula to rotate a zonal harmonic to any direction.

The forumal can be found in http://www.cs.dartmouth.edu/~wjarosz/publications/dissertation/appendixB.pdf B.37.

The constants can be deducted from all of this. It involves the SH representation of a cosinus function around z axis, and some term depending on that comes out of B.37.

[EDIT] I forgot to mention, you're actually not just computing SH coefficients for light but you're actually integrating a light source with a material's BRDF so you're convolving their SH representation too. You need to take that into account when deducting constants.

Edited by vlj

##### Share on other sites

Okay. So I have come up with this for the lighting function:

vector<vec3> ProjectLightFunction_Dir(int bandsCount, const vec3& dir, const vec3& col)
{
vector<float> shFuncs;

float r = sqrtf(dir.x*dir.x + dir.y*dir.y + dir.z*dir.z);
float theta = acosf(dir.z/r);
float phi = atanf(dir.y/dir.x);

for (int l = 0; l < bandsCount; l++)
{
for (int m = -l; m <= l; m++)
{
shFuncs.push_back(SphericalHarmonic(l, m, theta, phi));
}
}

vector<vec3> coeffs;
coeffs.resize(bandsCount * bandsCount);

for (int i = 0; i < bandsCount * bandsCount; i++)
{
coeffs[i].x = 0.0f;
coeffs[i].y = 0.0f;
coeffs[i].z = 0.0f;
}

for (int i = 0; i < bandsCount * bandsCount; i++)
{
coeffs[i].x += col.x * (shFuncs[i]);
coeffs[i].y += col.y * (shFuncs[i]);
coeffs[i].z += col.z * (shFuncs[i]);
}

return coeffs;
}


It's just constant color function projected to SH in the direction of the light. Together with the "general" Oliveira code for Transfer function (pages 20-21) it works nicely.

I also integrated the cos(theta) function what yielded PI (https://www.wolframalpha.com/input/?i=int+%28+cos%28y%29*sin%28y%29+%29+dy+dx%2C+x%3D0+to+2*PI%2C+y%3D0+to+PI%2F2).

So I tried this transfer function code:

vector<vec3> ProjectTransferFunction_Dir(int bandsCount, const vec3& normal)
{
vector<float> shFuncs;

float r = sqrtf(normal.x*normal.x + normal.y*normal.y + normal.z*normal.z);
float theta = acosf(normal.z/r);
float phi = atanf(normal.y/normal.x);

for (int l = 0; l < bandsCount; l++)
{
for (int m = -l; m <= l; m++)
{
shFuncs.push_back(SphericalHarmonic(l, m, theta, phi));
}
}

vector<vec3> coeffs;
coeffs.resize(bandsCount * bandsCount);

for (int i = 0; i < bandsCount * bandsCount; i++)
{
coeffs[i].x = 0.0f;
coeffs[i].y = 0.0f;
coeffs[i].z = 0.0f;
}

for (int i = 0; i < bandsCount * bandsCount; i++)
{
coeffs[i].x += PI * (shFuncs[i]);
coeffs[i].y += PI * (shFuncs[i]);
coeffs[i].z += PI * (shFuncs[i]);
}

return coeffs;
}


But using this transfer function gives me incorrect results. The lighting is not only incorrect but also too intense.

##### Share on other sites

As I wrote before, Lambertian directional light is projected by evaluating the SH basis function in the light's direction. You can derive it by writing down SH projection integral and using orthogonality property (see Funk-Hecke theorem) to simplify. So at the end just take your light direction vector, plug into Ylm and scale results by light's color. Eg. for 3rd order SH:

shR[ 0 ] = +0.282095f * lightColor.r;
shR[ 1 ] = -0.488603f * n.y * lightColor.r;
shR[ 2 ] = +0.488603f * n.z * lightColor.r;
shR[ 3 ] = -0.488603f * n.x * lightColor.r;
...
shG[ 0 ] = +0.282095f * lightColor.g;
...
shB[ 0 ] = +0.282095f * lightColor.b;
...

Edited by Krzysztof Narkowicz

##### Share on other sites

Is n in the formula the light direction vector? If so, I understand those formulas you wrote down. This is simply projection of lightColor function to SH in the light direction. This is also what the code function ProjectLightFunction_Dir does.

What bugs me is why the final lighting is computed like this:

     vec3 norm = gl_Normal;

color  = Coefficients[0];
color += Coefficients[1] * norm.x;
color += Coefficients[2] * norm.y;
color += Coefficients[3] * norm.z;
color += Coefficients[4] * norm.x * norm.z;
color += Coefficients[5] * norm.y * norm.z;
color += Coefficients[6] * norm.x * norm.y;
color += Coefficients[7] * (3.0 * norm.z * norm.z  - 1.0);
color += Coefficients[8] * (norm.x * norm.x  - norm.y * norm.y);


Why are coefficients of the transfer function (?) projected to SH the way they are? How are they computed?

Given that the code above is correct it appears that, for instance, coefficient c of the transfer function is: c(l=2, m=-1) = y*z. But why is it so?

Edited by maxest

##### Share on other sites

One more question.

I found this code:

for each texel of the cubemap
{
float3 c_light = texel_radiance;
float weight = texelSolidAngle;
float3 n = texelDirection;
float SHLightL[9];
SHLightL[0] = 0.282095f * c_light * weight;
SHLightL[1] = -0.488603f * n.y * c_light * weight;
SHLightL[2] = 0.488603f * n.z * c_light * weight;
SHLightL[3] = -0.488603f * n.x * c_light * weight;
SHLightL[4] = 1.092548f * n.x * n.y * c_light * weight;
SHLightL[5] = -1.092548f * n.y * n.z * c_light * weight;
SHLightL[6] = 0.315392f * (3.0f * n.z * n.z - 1.0f) * c_light * weight;
SHLightL[7] = -1.092548f * n.x * n.z * c_light * weight;
SHLightL[8] = 0.546274f * (n.x * n.x - n.y * n.y) * c_light * weight;
}
// Convolve with SH-projected cosinus lobe
float ConvolveCosineLobeBandFactor[] =
{
PI,
2.0f * PI/3.0f, 2.0f * PI/3.0f, 2.0f * PI/3.0f,
PI/4.0f, PI/4.0f, PI/4.0f, PI/4.0f, PI/4.0f
}
for (int i = 0; i < 9; ++i)
SHLightL[i] *= ConvolveCosineLobeBandFactor[i];


Now how are ConvolveCosineLobeBandFactor calculated?

To get the first one, I tried solving integral of ( cos(theta) * sin(theta) * Y(l=0, m=0, theta, phi) dtheta dphi ) but didn't end up with PI.