I am currently trying to verify the correctness of my diffuse lighting. I am fairly sure that my GGX specular model is correct, but I'm struggling with diffuse ambient lighting.

I am using a normalized lambert diffuse model. For simplicity, I'm assuming that the light has a color/intensity of (1, 1, 1) and that the surface has an albedo of (1, 1, 1) as those can be trivially multiplied into the final result.

float diffuse = saturate(dot(N, L)) / PI;

QUESTION 1: This implies that the correct result for a directional light hitting the surface at a perfect 90 degree angle (meaning the dot product is 1.0), I'd get a result of 1/PI. This is indeed what I'm seeing. Is that the expected/correct result?

Secondly, I'm looking into my usage of irradiance maps. I currently generate a tiny irradiance cube map, which I do a cosine convolution of. I believe the result of this convolution is correct. I am currently simply sampling this map and applying that as a diffuse color. For simplicity's sake, let's assume that the environment map is (1, 1, 1) in all directions (both before and after convolution).

vec3 ambient = texture(irradianceMap, N);

QUESTION 2: What is the expected amount of diffuse light reflected off of a white object in this case? In my current code, I simply sample the irradiance map and get 1.0, but online sources are saying that the result in this case should be 1/PI as well.

For various reasons, I want to move over to storing my irradiance map with spherical harmonics instead. I have a lot of questions about the math of doing ambient lighting with SH, specifically about how the math just all seems to cancel out, but I'd like to know what the expected diffuse light should be before I proceed with that.

I'm haunted by the same uncertainty about this damn division by PI, since many years.
I assume it comes from this: If we apply cosine weighting to a half sphere, representing the incoming light, it ‘flattens’ the sphere and it becomes a disk, so we need to divide by it's area of PI to ‘normalize’ the integrated, incoming light sampled over the hemisphere.

However, such speculation does not help me against the uncertainty. So i always need to trial and error and compare with reference images made with a pathtracer.

But i can tell there is no need to wait with your SH work. Whatever you do, the division by PI is usually just a constant factor affecting the end result, to me at least. So you can decide if it's needed or not later. It should not really affect any work.

Thanks for the response!

I believe the normalization factor for lambertian diffuse comes from doing a double integral over all possible directions in the relevant hemisphere (where dot(N, L) is positive) yields 0.5, while the normalized average should be 1/(2pi), one over the area of a hemisphere. Therefore, the normalization factor becomes (what we want) / (what we get) = (1/(2pi)) / 0.5 = 1/pi. I did the math at some point and I believe it to be correct. However, I'm not sure if the same math applies when sampling an irradiance map.

The main reason I want to make sure my irradiance maps work as they should before throwing them out is to have something to compare with when I implement SH irradiance maps.

When investigating SH ambience, I've been trying to simplify the math and make it “make sense”. Looking again at an irradiance map of intensity (1, 1, 1), all but the first term will be 0, so a 0th order "SH" would be enough in this case.

Following the math, the order 0 coefficient is ½*sqrt(1/pi).

From online sources, it looks like this is both the coefficient when sampling the radiance map and writing out the coefficients AND the coefficient when evaluating the SH in a certain direction. In addition, there's a second coefficient when sampling, A0, which is equal to pi.

So the effective coefficient is: pi * (0.5*sqrt(1/pi))^2 = 0.25. We'd then possibly need to divide by pi again here.

The problem is that 0.25 is “wildly” different from the expected value of 1.0, as I would expect the same result as if I evaluated the irradiance map directly. Somehow the math doesn't seem to work out to what I expect here. This is what made me partly doubt my original irradiance map code in the first place

For the 1st order SH coefficients, we get that the coefficient is 0.5*sqrt(3/pi), and A1 = 2pi/3. Simplifying (0.5*sqrt(3/pi))^2 * 2pi/3 yields 0.5. Through experiments, I found that an effective coefficient of 2.0 yields the correct result, again off by exactly a factor of 4.

What's going on with that? Why are my SH calculations off by a factor of 4?

theagentd said:
I did the math at some point and I believe it to be correct.

Sounds very similar to what i have meant.

theagentd said:
I'm not sure if the same math applies when sampling an irradiance map.

It depends on what exactly is in your map. Maybe you already did the division on the data in the map, maybe you didn't and need to do it after sampling from it.
A matter of conventions and details, thus it's probably hard to clarify. We are not the only ones struggling with this.

theagentd said:
The main reason I want to make sure my irradiance maps work as they should before throwing them out is to have something to compare with when I implement SH irradiance maps.

Well, my story is quite similar: I work on realtime GI, and at some point i wanted to compare with ground truth. So i worked on path tracer for the first time. It was easy. I could implement it without looking up any tutorials on path tracing. There was only one uncertainty: Do i need to divide by PI or not? I've asked her on the forum about it, and got some help, but i could not be sure about misconceptions.
So, instead using the path tracer to validate my realtime GI, i did it the other way around. Doing the division the results were almost equivalent, so i took that as a yes for the question, which also agreed with the answers i got here. Later i could confirm it 100% by comparing with a commercial offline rendering tool.

theagentd said:
The problem is that 0.25 is “wildly” different from the expected value of 1.0

I may be wrong, but i could swear the same confusion happened to me as well when working with SH. :D
Sadly i don't have the code around currently, but i may remember a multiplication by 4 which i felt pretty confused about, but it worked as expected with it. (not sure if the same factor worked for both SH2 and SH3)

But again i can't help. I'm completely self thought with math, thus happy if i can use such things at all, without understanding them fully.

Do you know any guy who really understands how quaternion multiplication works?
I don't. Really nobody does. And still we all use it without worries, because it just works.

You may get better responses, but if not - just multiply by 4 and call it a day. Maybe one day you'll figure out the exact reason of why it's needed… \:D/

theagentd said:
What's going on with that? Why are my SH calculations off by a factor of 4?

The basic process for using SH is to first have some data in spatial domain (i.e. cube map), then project it into SH domain. For 0-order SH that simply involves computing a weighted average of the input values over all directions, where the weights are all 0.5*sqrt(1/pi). There is an additional normalization factor of 4pi that is also applied at this stage, which is the result of integrating the directions over the unit sphere.

Then, to use the SH values you first evaluate the SH function for the direction you are interested in, then compute the dot product of those values with the SH coefficients produced by the projection. For 0-order SH, this means multiplying again by 0.5*sqrt(1/pi), resulting in 1/(4*pi). This 1/(4pi) cancels out the 4pi that we multiplied with the SH data for normalization, so we get out the same values as the input cube map.

Here is code illustrating the process. This is derived from code that I know is correct.

H; // space domain input data (eg. cube map)
H_lm = 0; // accumulator containing SH data

for ( int i = 0; i < directionCount; i++ )
{
    // Evaluate the spherical harmonic basis functions for the sample direction.
    Y_lm = evaluateSH( order, directions[i] );
    
    // Multiply input data for direction with SH for direction, and accumulate.
    H_lm += Y_lm .* H( directions[i] );
}

// Normalize based on the number of samples and unit sphere surface area.
const float normalize = (4.0f*PI) / directionCount;
H_lm *= normalize;

//-------------------------

// Interpolate SH coefficients
Y_lm = evaluateSH( order, my_random_direction );
my_random_direction_data = sum( Y_lm .* H_lm );

Note that this is pseudocode. The “.*” operator is an element-wise multiplication of an array.

Aressera said:
There is an additional normalization factor of 4pi that is also applied at this stage, which is the result of integrating the directions over the unit sphere.

That's what i remembered.

But do you eventually know if there is a way to find local minima in higher order SH, using an analytical method?

I worked on this, but could not find a way. So i ended up with iterative search, but that's not so nice.

@joej

I've studied a bit of math at uni, but not this kind of math. I've got a pretty good intuition of frequency domain transformations, so SHs make intuitive sense to me as a Fourier transform over the surface of a sphere. It's just that the math is so complex and unintuitive that I'm losing my miiiiind.

I appreciate the help! I'm just hoping to build some more intuition about where the 4 is coming from. 🙂

Aressera said:

Here is code illustrating the process. This is derived from code that I know is correct.

Interesting. I agree with the intuition that the 0-order SH should just be a weighted average of the entire sphere. That pseudocode makes a lot of sense and my code is quite similar. The only thing different is the 4pi normalization factor instead of just pi.

From what I can tell, there's a different normalization factor for each order:

• A0 = pi
• A1 = 2pi/3
• A2 = pi/4
• A4 = -pi/24
• …

It looks like you're using a factor of 4pi for all of them. May I ask where this factor comes from? It'd also be very useful to see how your evaluateSH() function works, given that some of the factors may be multiplied in there.

Aressera said:

For 0-order SH, this means multiplying again by 0.5*sqrt(1/pi), resulting in 1/(4*pi). This 1/(4pi) cancels out the 4pi that we multiplied with the SH data for normalization, so we get out the same values as the input cube map.

Why the hell are we even doing these weird calculations in the first place, if it's all gonna cancel out in the end anyway? We should be able to bake in all constants from the normalization factor(s) and evaluateSH(), leaving only the parts that are a function of x, y and z in there. That is what I was trying to do with my math, but kept ending up off by a factor of 4. Why are people even doing it this clearly inefficient way? You'd want to bake as much into the precomputed SHs as possible, so that evaluation is cheap.

Final question for @aressera: Do you apply the 1/pi factor when applying the resulting irradiance as ambient light?

theagentd said:
From what I can tell, there's a different normalization factor for each order:

The constants are different, but the normalization is the same (1/sqrt(4pi) in all cases). The difference is that in higher order SH the values are split over more and more coefficients so the constants are different. If you compute the inner product of SH coefficients you get 1/(4pi) out, for any order:

const Size order = 3;
const Size count = (order + 1)*(order + 1);

PODArray<Float32> sh;
PODArray<Float32> sh2;
sh.allocate( count );
sh2.allocate( count );

RandomFloat32 rand;

Float64 total = 0.0f; // double required
const Size N = 100000000;
for ( Index j = 0; j < N; j++ )
{
    Vector3f v = sampleSphereUniform( rand.sample11<Float32>(), rand.sample01<Float32>() );
    v = v.normalize();
    math::SH::evaluate( order, v.x, v.y, v.z, sh.getPointer() + 0 );
    
    Vector3f v2 = sampleSphereUniform( rand.sample11<Float32>(), rand.sample01<Float32>() );
    v2 = v2.normalize();
    math::SH::evaluate( order, v2.x, v2.y, v2.z, sh2.getPointer() + 0 );
    
    for ( Index i = 0; i < count; i++ )
        total += sh[i]*sh2[i];
}

// This always prints something close to 1/(4pi), for any order
Console << total / N;

As I said in the original post, the 1/(4pi) comes from the integral over a unit sphere, and is the same as 1/(surface area) for a unit sphere.

I use the fast SH evaluation functions from Peter-Pike Sloan, see paper and code “Efficient Spherical Harmonic Evaluation” (2013). Here is the generated function for order 2:

template < typename T, typename U >
void evaluateSH( U x, U y, U z, T* expansion )
{
	U fX = x, fY = y, fZ = z;
	U fC0, fS0, fTmpA;
	U fC1, fS1, fTmpB, fTmpC;
	U fZ2 = fZ*fZ;
	
	expansion[0] = T(U(0.2820947917738781));
	expansion[2] = T(U(0.4886025119029199)*fZ);
	expansion[6] = T(U(0.9461746957575601)*fZ2 + U(-0.31539156525252));
	fC0 = fX;
	fS0 = fY;
	
	fTmpA = U(0.48860251190292);
	expansion[3] = T(fTmpA*fC0);
	expansion[1] = T(fTmpA*fS0);
	fTmpB = (U(1.092548430592079)*fZ);
	expansion[7] = T(fTmpB*fC0);
	expansion[5] = T(fTmpB*fS0);
	fC1 = fX*fC0 - fY*fS0;
	fS1 = fX*fS0 + fY*fC0;
	
	fTmpC = U(0.5462742152960395);
	expansion[8] = T(fTmpC*fC1);
	expansion[4] = T(fTmpC*fS1);
}

It won't get much more efficient than that.

theagentd said:
Why the hell are we even doing these weird calculations in the first place, if it's all gonna cancel out in the end anyway?

Because we need to ensure the SH are normalized so that they behave correctly, otherwise you would get out the wrong value when interpolating the SH data, or some other property of SH might not work (e.g. rotation). You can apply the normalization at any step (e.g. during projection like in my pseudocode), but it needs to be there for mathematical correctness. Multiplying by 4pi is not “inefficient” as you say, there are no extra calculations in your shader.

theagentd said:
Final question for @aressera: Do you apply the 1/pi factor when applying the resulting irradiance as ambient light?

I checked my code and I'm not, but that is an error for historical reasons and I plan to fix it. You need to normalize the BRDF by 1/pi for a similar reason the SH are normalized to 1/sqrt(4pi). You can treat the irradiance cube map as just another directional light source, and use the same BRDF calculations as other light sources.

I've confirmed that your evaluateSH() function produces the same results as my test code (barring floating point errors):

		coefficients[0] = 0.5 * sqrt(1/PI);
		
		coefficients[1] = 0.5 * sqrt(3/PI) * y;
		coefficients[2] = 0.5 * sqrt(3/PI) * z;
		coefficients[3] = 0.5 * sqrt(3/PI) * x;

		coefficients[4] = 0.5 * sqrt(15/PI) * x*y;
		coefficients[5] = 0.5 * sqrt(15/PI) * y*z;
		coefficients[6] = 0.25 * sqrt(5/PI) * (3.0*z*z - 1.0);
		coefficients[7] = 0.5 * sqrt(15/PI) * x*z;
		coefficients[8] = 0.25 * sqrt(15/PI) * (x*x - y*y);

However, I need to use a different normalization factor than you. I need to use A0=4pi, A1=8pi/3 and A2=pi, which is 4 times the constants listed here: https://patapom.com/blog/SHPortal/#using-the-result-for-real-time-diffuse-irradiance-estimation​

		float a0 = 4*PI;
		float a1 = 8*PI/3;
		float a2 = PI;
		
		normalization[0] = a0;
		
		normalization[1] = a1;
		normalization[2] = a1;
		normalization[3] = a1;

		normalization[4] = a2;
		normalization[5] = a2;
		normalization[6] = a2;
		normalization[7] = a2;
		normalization[8] = a2;

I have a small brute-force Shadertoy shader for verifying this. Left half shows brute-force irradiance calculation, right half shows the same using 2nd order SHs, Here's the result using the a0, a1 and a2 variables:

If I instead use PI*4 for everything, I get the following:

Unfortunately, I don't see how a 4pi normalization factor can be correct for all orders. :(

Aressera said:

theagentd said:
Why the hell are we even doing these weird calculations in the first place, if it's all gonna cancel out in the end anyway?

Because we need to ensure the SH are normalized so that they behave correctly, otherwise you would get out the wrong value when interpolating the SH data, or some other property of SH might not work (e.g. rotation). You can apply the normalization at any step (e.g. during projection like in my pseudocode), but it needs to be there for mathematical correctness. Multiplying by 4pi is not “inefficient” as you say, there are no extra calculations in your shader.

Looking at the code, we basically have a sum over a large number of directions. For the 0th order SH, we basically have this:

sum( 0.5 * sqrt(1/PI) * colors );

We can trivially move the constants outside of the sum

( 0.5 * sqrt(1/PI) ) * sum(colors)

When evaluating the SH, we then multiply this by the same constant again, and by the normalization factor:

( 0.5 * sqrt(1/PI) ) * ( 0.5 * sqrt(1/PI) ) * ( 4.0 * PI ) * sum(colors)

This all simplifies to just:

sum(colors)

Similarly, all the 1st order SHs simplify to this (assuming A1=8pi/3 used as normalization factor):

2.0 * sum(colors * x|y|z)

It seems like we're doing a bunch of redundant calculations here that all just cancel out. These simple constants would not break any interpolation, as a simple scaling factor can be applied both before or after the interpolation. I can't really see how just storing these much simpler sums would break anything. :/

theagentd said:
I can't really see how just storing these much simpler sums would break anything

You might be misunderstanding something if you think anything can be precomputed beyond the SH projection H_lm from my previous pseudocode example. Interpolating the SH projection is as simple as evaluating the SH for the direction you are interested in, and computing the dot product with the SH data. You can't precompute anything there unless you restrict yourself to only specific directions.

If you use alternative normalization schemes (e.g. as used in ambisonics audio, commonly to fit better into fixed point), then you have to compensate your SH decoding function in the shader with the inverse of the normalization factors, so that when you calculate the dot product you get the same result as before without the normalization.