The following is taken from my current understanding of spherical harmonics. There may be errors involved which solve my question if they're resolved, so please notify me if I misunderstood something.
Evaluating reflected diffuse indirect illumination at a surface point is done by computing the integral (simplified rendering equation) in the last line of
(stolen from here).
An important property of functions in SH representation is that the integral in the last line boils down to a dot product between the vectors containing the SH coefficients for projected functions L(x,omega_i) and max(N_x dot omega_i,0).
For example, for each surface point in our target image that we wish to compute indirect illumination on, L(x,omega_i) is available in a precomputed vector of a 4 elements (for 2 bands SH approximation) containing spherical harmonics coefficients. It's not really important where it came from, just that it's available at this point.
From my understanding, to actually evaluate the lighting at a surface point [x] is to project the function max(N_x dot omega_i,0) into spherical harmonics, and dot it with the coefficient vector of the incoming radiance (and then, multiply it by the stuff in front of the integral, surface albedo / pi). The first 4 coefficients (2 bands) for max(N_x dot omega_i,0) projected into SH basis are known, they are:
c_0 = 0.88622692545...
c_1 = -1.02332670795 * n_y
c_2 = 1.02332670795 * n_z
c_3 = -1.02332670795 * n_x
So the irradiance at that surface point with normal n would just be
E = dot(L_sh,n_sh)
where L_sh is the coefficient vector for the incoming radiance, and n_sh the coefficient vector for the cosine lobe that we just calculated. Specifically, L_sh is the same max(...) function projected into SH basis [with the normal from which the light originated], just scaled a little by the color its transporting.
The above actually looks fine and seems to produce somewhat correct lighting (hard to judge without ground truth reference). What's confusing me now, is that I've come across several sources (for example here, page 17), which, in the final step, when evaluating the lighting, don't actually use the coefficients from above to project max(N_x dot omega_i,0) into SH basis, but instead calculate n_sh as
n_sh0 = 0.282094792f
n_sh1 = -0.488602512f * -n_y
n_sh2 = 0.488602512f * -n_z
n_sh3 = -0.488602512f * -n_x
These four coefficients are actually just the spherical harmonics basis function coefficients parametrized in cartesian coordinates (here they are for spherical coordinates).
What's up with that? Why aren't they projecting the max(...) function into SH basis instead? Am I misunderstanding something? Unfortunately, the usual sources (gritty details and stupid SH tricks) don't seem to go into details of this either.