Analytical projection of various light types into L2 Spherical Harmonics. https://pema.dev/obsidian/math/light-transport/analytically-projecting-lights-into-sh.html

AnalyticalSHProjection.cginc

// SPDX-License-Identifier: MIT
// Author: pema99

// This file contains various functions for projecting irradiance due to contribution from 
// various light types into L2 spherical harmonics.

// Quick reference:
// Functions for projecting irradiance:
//   - SphericalHarmonicsL2 ProjectDirectionalLightIrradianceSHL2(lightDirection, lightColor)
//   - SphericalHarmonicsL2 ProjectPointLightIrradianceSHL2(shadingPosition, lightPosition, lightColor)
//   - SphericalHarmonicsL2 ProjectSpotLightIrradianceSHL2(shadingPosition, lightPosition, lightDirection, lightColor, outerHalfAngle, innerHalfAngle)
//   - SphericalHarmonicsL2 ProjectSphereLightIrradianceSHL2(shadingPosition, lightPosition, lightRadius, lightColor)
//   - SphericalHarmonicsL2 ProjectQuadLightIrradianceSHL2(shadingPosition, lightVertices, lightColor)
// Helpers functions for evaluating SH:
//   - float3 EvaluateSHL2(sh, direction)
//   - void SHL2ToUnityUniforms(sh, out SHAr, out SHAg, out SHAb, out SHBr, out SHBg, out SHBb, out SHC)

namespace ASH
{
    struct SphericalHarmonicsL2
    {
        float3 coefficients[9];
    };

    #ifndef PI
    static const float PI = 3.14159265358979323846;
    #endif

    // Constant parts of each SH basis functions
    static const float L0_NORMALIZATION = 0.5 * sqrt(1.0 / PI);
    static const float L1_NORMALIZATION = 0.5 * sqrt(3.0 / PI);
    static const float L2_2_NORMALIZATION = 0.5 * sqrt(15.0 / PI);
    static const float L20_NORMALIZATION = 0.25 * sqrt(5.0 / PI);
    static const float L22_NORMALIZATION = 0.25 * sqrt(15.0 / PI);

    // Constants for radiance -> irradiance convolution
    static const float A_HAT_0 = PI;
    static const float A_HAT_1 = 2.0 * PI / 3.0;
    static const float A_HAT_2 = PI / 4.0;

    // Number of vertices for a polygon light.
    static const uint NUM_VERTS = 4;

    // SH basis functions
    float Y0() // Basis function l=0
    {
        return L0_NORMALIZATION;
    }

    float Y1_1(float3 dir) // Basis function l=1, m=-1
    {
        return L1_NORMALIZATION * dir.y;
    }

    float Y10(float3 dir) // Basis function l=1, m=0
    {
        return L1_NORMALIZATION * dir.z;
    }

    float Y11(float3 dir) // Basis function l=1, m=1
    {
        return L1_NORMALIZATION * dir.x;
    }

    float Y2_2(float3 dir) // Basis function l=2, m=-2
    {
        return L2_2_NORMALIZATION * dir.x * dir.y;
    }

    float Y2_1(float3 dir) // Basis function l=2, m=-1
    {
        return L2_2_NORMALIZATION * dir.y * dir.z;
    }

    float Y20(float3 dir) // Basis function l=2, m=0
    {
        return L20_NORMALIZATION * (3.0 * dir.z * dir.z - 1);
    }

    float Y21(float3 dir) // Basis function l=2, m=1
    {
        return L2_2_NORMALIZATION * dir.x * dir.z;
    }

    float Y22(float3 dir) // Basis function l=2, m=2
    {
        return L22_NORMALIZATION * (dir.x * dir.x - dir.y * dir.y);
    }

    // Associated legendre polynomial for zonal harmonics, up to L3
    float AssociatedLegendreZH(uint l, float x)
    {
        [forcecase] switch (l)
        {
            case 0: return 1.0;
            case 1: return x;
            case 2: return 0.5 * (3.0 * x * x - 1.0);
            case 3: return 0.5 * x * (5.0 * x * x - 3.0); 
            default: return 0.0;
        }
    }

    // Normalization factor for zonal harmonics
    float NormalizationZH(uint l)
    {
        return sqrt((2.0 * l + 1.0) / (4.0 * PI));
    }

    // Evaluate L2 spherical harmonics at a given direction
    float3 EvaluateSHL2(SphericalHarmonicsL2 sh, float3 direction)
    {
        float3 sum = 0.0;

        sum += sh.coefficients[0] * Y0();

        sum += sh.coefficients[1] * Y1_1(direction);
        sum += sh.coefficients[2] * Y10(direction);
        sum += sh.coefficients[3] * Y11(direction);
        
        sum += sh.coefficients[4] * Y2_2(direction);
        sum += sh.coefficients[5] * Y2_1(direction);
        sum += sh.coefficients[6] * Y20(direction);
        sum += sh.coefficients[7] * Y21(direction);
        sum += sh.coefficients[8] * Y22(direction);

        return sum;
    }

    // Convert SHL2 coefficients to Unity's uniform format
    void SHL2ToUnityUniforms(
        SphericalHarmonicsL2 sh,
        out float4 SHAr,
        out float4 SHAg,
        out float4 SHAb,
        out float4 SHBr,
        out float4 SHBg,
        out float4 SHBb,
        out float4 SHC)
    {
        // Multiply constant part of each basis function
        sh.coefficients[0] *= L0_NORMALIZATION;
        sh.coefficients[1] *= L1_NORMALIZATION;
        sh.coefficients[2] *= L1_NORMALIZATION;
        sh.coefficients[3] *= L1_NORMALIZATION;
        
        sh.coefficients[4] *= L2_2_NORMALIZATION;
        sh.coefficients[5] *= L2_2_NORMALIZATION;
        sh.coefficients[6] *= L20_NORMALIZATION;
        sh.coefficients[7] *= L2_2_NORMALIZATION;
        sh.coefficients[8] *= L22_NORMALIZATION;

        // Map to shader uniforms
        SHAr = float4(sh.coefficients[3].r, sh.coefficients[1].r, sh.coefficients[2].r, sh.coefficients[0].r - sh.coefficients[6].r);
        SHAg = float4(sh.coefficients[3].g, sh.coefficients[1].g, sh.coefficients[2].g, sh.coefficients[0].g - sh.coefficients[6].g);
        SHAb = float4(sh.coefficients[3].b, sh.coefficients[1].b, sh.coefficients[2].b, sh.coefficients[0].b - sh.coefficients[6].b);

        SHBr = float4(sh.coefficients[4].r, sh.coefficients[5].r, sh.coefficients[6].r * 3.0f, sh.coefficients[7].r);
        SHBg = float4(sh.coefficients[4].g, sh.coefficients[5].g, sh.coefficients[6].g * 3.0f, sh.coefficients[7].g);
        SHBb = float4(sh.coefficients[4].b, sh.coefficients[5].b, sh.coefficients[6].b * 3.0f, sh.coefficients[7].b);
    
        SHC = float4(sh.coefficients[8].r, sh.coefficients[8].g, sh.coefficients[8].b, 1.0);
    }

    // Project irradiance from directional light into L2 spherical harmonics
    SphericalHarmonicsL2 ProjectDirectionalLightIrradianceSHL2(
        float3 lightDirection,
        float3 lightColor)
    {
        SphericalHarmonicsL2 sh;
    
        float3 toLightNorm = normalize(-lightDirection);
                    
        sh.coefficients[0] = lightColor * Y0() * A_HAT_0;
    
        sh.coefficients[1] = lightColor * Y1_1(toLightNorm) * A_HAT_1;
        sh.coefficients[2] = lightColor * Y10(toLightNorm)  * A_HAT_1;
        sh.coefficients[3] = lightColor * Y11(toLightNorm)  * A_HAT_1;
    
        sh.coefficients[4] = lightColor * Y2_2(toLightNorm) * A_HAT_2;
        sh.coefficients[5] = lightColor * Y2_1(toLightNorm) * A_HAT_2;
        sh.coefficients[6] = lightColor * Y20(toLightNorm)  * A_HAT_2;
        sh.coefficients[7] = lightColor * Y21(toLightNorm)  * A_HAT_2;
        sh.coefficients[8] = lightColor * Y22(toLightNorm)  * A_HAT_2;

        return sh;
    }

    // Project irradiance from point light into L2 spherical harmonics
    SphericalHarmonicsL2 ProjectPointLightIrradianceSHL2(
        float3 shadingPosition,
        float3 lightPosition,
        float3 lightColor)
    {
        SphericalHarmonicsL2 sh;
    
        float3 toLight = lightPosition - shadingPosition;
        float distanceSquared = dot(toLight, toLight);
        float3 toLightNorm = normalize(toLight);
                    
        sh.coefficients[0] = lightColor / distanceSquared * Y0() * A_HAT_0;
    
        sh.coefficients[1] = lightColor / distanceSquared * Y1_1(toLightNorm) * A_HAT_1;
        sh.coefficients[2] = lightColor / distanceSquared * Y10(toLightNorm)  * A_HAT_1;
        sh.coefficients[3] = lightColor / distanceSquared * Y11(toLightNorm)  * A_HAT_1;
    
        sh.coefficients[4] = lightColor / distanceSquared * Y2_2(toLightNorm) * A_HAT_2;
        sh.coefficients[5] = lightColor / distanceSquared * Y2_1(toLightNorm) * A_HAT_2;
        sh.coefficients[6] = lightColor / distanceSquared * Y20(toLightNorm)  * A_HAT_2;
        sh.coefficients[7] = lightColor / distanceSquared * Y21(toLightNorm)  * A_HAT_2;
        sh.coefficients[8] = lightColor / distanceSquared * Y22(toLightNorm)  * A_HAT_2;

        return sh;
    }

    // Project irradiance from spot light into L2 spherical harmonics
    SphericalHarmonicsL2 ProjectSpotLightIrradianceSHL2(
        float3 shadingPosition,
        float3 lightPosition,
        float3 lightDirection,
        float3 lightColor,
        float outerHalfAngle,
        float innerHalfAngle)
    {
        SphericalHarmonicsL2 sh;
    
        float3 toLight = lightPosition - shadingPosition;
        float distanceSquared = dot(toLight, toLight);
        float3 toLightNorm = normalize(toLight);

        float cosTheta = dot(toLightNorm, -lightDirection);
        float theta = acos(clamp(cosTheta, -1, 1));
        float falloff = smoothstep(outerHalfAngle, innerHalfAngle, theta);

        sh.coefficients[0] = lightColor / distanceSquared * falloff * Y0() * A_HAT_0;
    
        sh.coefficients[1] = lightColor / distanceSquared * falloff * Y1_1(toLightNorm) * A_HAT_1;
        sh.coefficients[2] = lightColor / distanceSquared * falloff * Y10(toLightNorm)  * A_HAT_1;
        sh.coefficients[3] = lightColor / distanceSquared * falloff * Y11(toLightNorm)  * A_HAT_1;
    
        sh.coefficients[4] = lightColor / distanceSquared * falloff * Y2_2(toLightNorm) * A_HAT_2;
        sh.coefficients[5] = lightColor / distanceSquared * falloff * Y2_1(toLightNorm) * A_HAT_2;
        sh.coefficients[6] = lightColor / distanceSquared * falloff * Y20(toLightNorm)  * A_HAT_2;
        sh.coefficients[7] = lightColor / distanceSquared * falloff * Y21(toLightNorm)  * A_HAT_2;
        sh.coefficients[8] = lightColor / distanceSquared * falloff * Y22(toLightNorm)  * A_HAT_2;

        return sh;
    }

    // Project irradiance from spherical light into L2 spherical harmonics
    // NOTE: Unity's emissives are too dark by a factor of PI. To match them, divide the output coefficients by PI.
    SphericalHarmonicsL2 ProjectSphereLightIrradianceSHL2(
        float3 shadingPosition,
        float3 lightPosition,
        float lightRadius,
        float3 lightColor)
    {
        SphericalHarmonicsL2 sh;
    
        float3 toLight = lightPosition - shadingPosition;
        float distanceSquared = dot(toLight, toLight);
        float3 toLightNorm = normalize(toLight);
    
        // sin = opposite / hypotenuse = radius / distance
        float sinAngle = lightRadius / max(sqrt(distanceSquared), 0.000001f);
        // sin^2 + cos^2 = 1  =>  cos^2 = 1 - sin^2  =>  cos = sqrt(1 - sin^2)
        float cosAngle = sqrt(saturate(1.0f - sinAngle * sinAngle)); 

        // Calculate projection in a space where +Z is the direction towards the light
        float l0Hat = sqrt(PI) * (1.0 - cosAngle);
        float l1Hat = sqrt(PI / 3.0) * (AssociatedLegendreZH(0, cosAngle) - AssociatedLegendreZH(2, cosAngle));
        float l2Hat = sqrt(PI / 5.0) * (AssociatedLegendreZH(1, cosAngle) - AssociatedLegendreZH(3, cosAngle));
    
        // Rotate into global space, multiply on intensity to get final coefficients
        sh.coefficients[0] = lightColor * sqrt(4.0 * PI) * Y0() * l0Hat * A_HAT_0;
    
        sh.coefficients[1] = lightColor * sqrt(4.0 * PI / 3.0) * Y1_1(toLightNorm) * l1Hat * A_HAT_1;
        sh.coefficients[2] = lightColor * sqrt(4.0 * PI / 3.0) * Y10(toLightNorm)  * l1Hat * A_HAT_1;
        sh.coefficients[3] = lightColor * sqrt(4.0 * PI / 3.0) * Y11(toLightNorm)  * l1Hat * A_HAT_1;

        sh.coefficients[4] = lightColor * sqrt(4.0 * PI / 5.0) * Y2_2(toLightNorm) * l2Hat * A_HAT_2;
        sh.coefficients[5] = lightColor * sqrt(4.0 * PI / 5.0) * Y2_1(toLightNorm) * l2Hat * A_HAT_2;
        sh.coefficients[6] = lightColor * sqrt(4.0 * PI / 5.0) * Y20(toLightNorm)  * l2Hat * A_HAT_2;
        sh.coefficients[7] = lightColor * sqrt(4.0 * PI / 5.0) * Y21(toLightNorm)  * l2Hat * A_HAT_2;
        sh.coefficients[8] = lightColor * sqrt(4.0 * PI / 5.0) * Y22(toLightNorm)  * l2Hat * A_HAT_2;
    
        return sh;
    }

    float3 BoundaryIntegral(float a, float b, float x)
    {
        float3 B_n;
        
        float z = a*cos(x) + b*sin(x);
        float tmp1 = a*sin(x) - b*cos(x);
        float tmp2 = a*a+b*b-1;
        
        B_n[0] = x;
        B_n[1] = tmp1 + b;

        float C_n = (tmp1 * z) + (tmp2 * x) + (B_n[0]) + (b * a);
        C_n *= 0.5;

        B_n[2] = 3.0 * C_n - B_n[0];
        B_n[2] *= 0.5;
        
        return B_n;
    }
                
    float3 EvalLight(float3 dir, float3 verts[NUM_VERTS], float3 gam[NUM_VERTS], float3 gamP[NUM_VERTS])
    {
        float3 total = 0;

        for (uint edge = 0; edge < NUM_VERTS; edge++)
        {
            uint next = (edge + 1) % 4;
            float3 bound = BoundaryIntegral(dot(dir, verts[edge]), dot(dir, gamP[edge]), acos(clamp(dot(verts[edge], verts[next]), -1, 1)));
            total += bound * dot(dir, gam[edge]);
        }
        
        return 0.5 * float3(0, total.x, total.y);
    }
                
    float SolidAngle(float3 verts[NUM_VERTS])
    {
        float sa = 0;
        for (uint i = 0; i < NUM_VERTS; i++)
        {
            uint next = (i + 1) % NUM_VERTS;
            uint prev = (i + NUM_VERTS - 1) % NUM_VERTS;
            float3 tmp1 = cross(verts[i], verts[prev]);
            float3 tmp2 = cross(verts[i], verts[next]);
            sa += acos(clamp(dot(tmp1, tmp2) / (length(tmp1) * length(tmp2)), -1, 1));
        }
        return sa - (NUM_VERTS - 2) * PI;
    }

    // Project irradiance from quad light (assumed to be planar) into L2 spherical harmonics
    // NOTE: Unity's area lights and emissives are too dark by a factor of PI. To match them, divide the output coefficients by PI.
    SphericalHarmonicsL2 ProjectQuadLightIrradianceSHL2(float3 shadingPosition, float3 lightVertices[NUM_VERTS], float3 lightColor)
    {
        SphericalHarmonicsL2 sh;
        for (uint i = 0; i < 9; i++)
        {
            sh.coefficients[i] = 0;
        }
        
        // Check which side of the light we are on (assumes planar polygon)
        float3 lightNormal = cross(lightVertices[2] - lightVertices[0], lightVertices[1] - lightVertices[0]);
        if (dot(lightNormal, shadingPosition - lightVertices[0]) < 0.0)
            return sh;
        
        for (uint i = 0; i < NUM_VERTS; i++)
        {
            lightVertices[i] = normalize(lightVertices[i] - shadingPosition);
        }
        
        float3 G[NUM_VERTS];
        float3 Gp[NUM_VERTS];
        for (int i = 0; i < NUM_VERTS; i++)
        {
            uint next = (i + 1) % NUM_VERTS;
            G[i] = normalize(cross(lightVertices[i], lightVertices[next]));
            Gp[i] = cross(G[i], lightVertices[i]);
        }

        float3 SA = SolidAngle(lightVertices);
        sh.coefficients[0] = lightColor * sqrt(1.0 / (4.0 * PI)) * SA;

        float3 w20 = EvalLight(float3(0.866025, -0.500001, -0.000004), lightVertices, G, Gp);
        float3 w21 = EvalLight(float3(-0.759553, 0.438522, -0.480394), lightVertices, G, Gp);
        float3 w22 = EvalLight(float3(-0.000002, 0.638694, 0.769461), lightVertices, G, Gp);
        float3 w23 = EvalLight(float3(-0.000004, -1.000000, -0.000004), lightVertices, G, Gp);
        float3 w24 = EvalLight(float3(-0.000007, 0.000003, -1.000000), lightVertices, G, Gp);

        sh.coefficients[1] = lightColor * dot(float3(2.1995339, 2.50785367, 1.56572711), float3(w20[1], w21[1], w22[1]));
        sh.coefficients[2] = lightColor * dot(float2(-1.82572523, -2.08165037), float2(w20[1], w21[1]));
        sh.coefficients[3] = lightColor * dot(float3(2.42459869, 1.44790525, 0.90397552), float3(w20[1], w21[1], w22[1]));
        
        sh.coefficients[4] = lightColor * dot(float3(-1.33331385, -0.66666684, -0.99999606), float3(w20[2], w23[2], w24[2]));
        sh.coefficients[5] = lightColor * dot(float3(1.1747938, -0.47923799, -0.69556433), float3(w22[2], w23[2], w24[2]));
        sh.coefficients[6] = lightColor * w24[2];
        sh.coefficients[7] = dot(float3(-1.21710396, 1.58226094, 0.67825711), float3(w20[2], w21[2], w22[2]));
        sh.coefficients[7] += dot(float2(-0.27666329, -0.76671491), float2(w23[2], w24[2]));
        sh.coefficients[7] *= lightColor;
        sh.coefficients[8] = lightColor * dot(float2(-1.15470843, -0.57735948), float2(w23[2], w24[2]));

        // Radiance to irradiance
        sh.coefficients[0] *= A_HAT_0;
        sh.coefficients[1] *= A_HAT_1;
        sh.coefficients[2] *= A_HAT_1;
        sh.coefficients[3] *= A_HAT_1;
        sh.coefficients[4] *= A_HAT_2;
        sh.coefficients[5] *= A_HAT_2;
        sh.coefficients[6] *= A_HAT_2;
        sh.coefficients[7] *= A_HAT_2;
        sh.coefficients[8] *= A_HAT_2;
        
        return sh;
    }
}