Sign in to follow this  
maxest

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 this post


Link to post
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 this post


Link to post
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 this post


Link to post
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 this post


Link to post
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 this post


Link to post
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 this post


Link to post
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 this post


Link to post
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 this post


Link to post
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.

Share this post


Link to post
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 this post


Link to post
Share on other sites

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

 

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 this post


Link to post
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 this post


Link to post
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 this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this