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?

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.

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?

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

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.

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

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?

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.