# Spherical Harmonics - analytical solution

## 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 [tt]dot(L, N) = L.x*N.x + L.y*N.y + L.z*N.z[/tt].

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
{
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.

##### Share on other sites
If you want to gain a better understanding of the cosine lobe convolution in SH, I would suggest reading this paper by Ravi Ramamoorthi and Pat Hanrahan. It covers the details of what Krzysztof and vlj touched on in their posts.

One thing to be aware of is that the code you just posted will compute the irradiance incident onto the surface, but not the diffuse reflectance. To compute Lambertian diffuse reflectance, you should multiply your irradiance by DiffuseAlbedo/Pi.

##### Share on other sites

Hm, after the two recent posts I understood, to some extent, what vlj talked about .

Looking at the paper you pointed to, MJP, I would assume that equation (19) should give me the result of (cos(theta) * Ylm) but this equation doesn't really yield PI, PI*2/3, PI/4.

I also tried integrating numerically the cosine lobe like this:

float CosineLobe0()
{
int n = 0;
float value = 0.0f;

for (float phi = 0.0f; phi <= 2.0f*PI; phi += 0.01f)
{
for (float theta = 0.0f; theta <= PI/2.0f; theta += 0.01f)
{
value += cosf(theta) * sinf(theta) * SphericalHarmonic(0, 0, theta, phi);
n++;
}
}

value /= n;

return value;
}


But I get... 0.089264. No matter what I think of I can't really turn that into PI.

Edited by maxest

##### Share on other sites

You should additionally multiply An by a normalization term (4*pi)/(2*l+1))^0.5. See eq. 23-24.

##### Share on other sites

Since the band index is 0:

0.089264 *  sqrt(4*pi) = 0.31643264109

Still not PI

##### Share on other sites

Eq. 19 has a small error. Instead of n!/2 it should be (n/2)!. So:

For even n: An = 2*pi*((2n + 1)/(4*pi))^0.5*(((-1)^(n/2 - 1))/((n + 2)*(n - 1)))*(n!/(2^n*((n/2)!)^2))
LambdaL = (4*pi)/(2*l+1))^0.5

A0 * Lambda0 ~= 0.886227 * 3.54491 = 3.14159

##### Share on other sites

Okay. Don't really understand the whole derivation of the formula, but the end result (equation (26)) seems to provide correct results.

Thanks :).

Could you also please give me a hint about what's wrong with my numerical integration of the cosine lobe? I suppose that simply dividing by n the final result is not correct. I probably should mul each "summation element" by the solid angle area but don't really know how to estimate it. Or maybe I'm doing something else wrong?

##### Share on other sites

I looked again at my numerical integration code I posted before and it's clear that it's severly flawed. Instead of integrating it's actually computing some sort of average value...

Anyway, I tried to integrate the function in Octave and it yielded 0.886227 indeed. I also fixed my code and here it is:

float CosineLobe0()
{
float step = 0.001f;
float value = 0.0f;

for (float phi = 0.0f; phi <= 2.0f*PI; phi += step)
{
float value2 = 0.0f;
for (float theta = 0.0f; theta <= PI/2.0f; theta += step)
{
float a = theta;
float b = theta + step;
float x = 0.5f * (a + b);

value2 += cosf(x) * sinf(x) * sphericalHarmonic(0, 0, theta, phi) * step;
}

value += value2 * step;
}

return value;
}


One last thing that bugs me now is equation (16). How come did spherical harmonic Y_n^0 become a one-argument function?

##### Share on other sites

Transfer function cos(phi) has no azimuthal dependency, so all terms for m != 0 will vanish after integration.

Thanks :).

## Create an account

Register a new account

• ### Forum Statistics

• Total Topics
628350
• Total Posts
2982211

• 10
• 9
• 24
• 11
• 9