I am building a 3d renderer in C++ using Vulkan.hpp. I have implemented Spherical Harmonics (SH) coefficients for diffuse radiance as part of global illumination. However, I am now trying to expand this functionality to adequately handle multiple probes. My SH calculation process and utilisation also uses a position cubemap as part of the depth testing while the SH generation uses an input environment map to generate the SH coefficients. To handle multiple probes, I have setup a cubemap array to handle the multiple position cubemaps for depth testing. These cubemaps could then theoretically be saved to the hard disk and only brought into the cubemap array when required. However, I'm not sure of the best way of handling SH coefficients for multiple probes. I want to ensure flexibility in the number of probes that are active in any given scene. What would be the best approach for this? To assist with context, I have provided the shader code to generate the SH coefficients and an excerpt of the lighting shader code where the coefficients are used:

SH_Ceofficient Compute Shader:

#version 450
#extension GL_EXT_scalar_block_layout : enable

#define PI 3.1415926535897932384626433832795
#define SH_ORDER 2
#define NUM_COEFFS ((SH_ORDER + 1) * (SH_ORDER + 1))
#define NUM_SAMPLES 64 // Or tune for quality/performance

layout(local_size_x = 1, local_size_y = 1, local_size_z = 1) in;

layout(scalar, set = 0, binding = 0) buffer SHBuffer {
    vec3 shCoeffs[];
} shBuffer;

layout(set = 0, binding = 1) uniform samplerCube irradianceMap;
layout(set = 0, binding = 2) uniform samplerCubeArray depthMap;

layout(set = 0, binding = 3) uniform ProbePosition {
    vec3 probePos;
} probePosition;

//layout(set = 0, binding = 4) uniform samplerCube normalMap;

// You may want to send other uniforms if needed, e.g. depth range parameters

// Real SH basis functions for 2nd order (9 coefficients), with same order as lighting shader
vec3 SHBasis(int i, vec3 dir) {
    float x = dir.x;
    float y = dir.y;
    float z = dir.z;

    // The SH basis functions closely follow the ones in the lighting shader:
    // Indexing same as lighting:
    // 0: 0.282095 (constant)
    // 1: 0.488603 * y
    // 2: 0.488603 * z
    // 3: 0.488603 * x
    // 4: 1.092548 * x * y
    // 5: 1.092548 * y * z
    // 6: 0.315392 * (3z^2 - 1)
    // 7: 1.092548 * x * z
    // 8: 0.546274 * (x^2 - y^2)

    // We return vec3 just to simplify multiply with radiance
    if (i == 0) return vec3(0.282095);
    else if (i == 1) return vec3(0.488603 * y);
    else if (i == 2) return vec3(0.488603 * z);
    else if (i == 3) return vec3(0.488603 * x);
    else if (i == 4) return vec3(1.092548 * x * y);
    else if (i == 5) return vec3(1.092548 * y * z);
    else if (i == 6) return vec3(0.315392 * (3.0 * z * z - 1.0));
    else if (i == 7) return vec3(1.092548 * x * z);
    else if (i == 8) return vec3(0.546274 * (x * x - y * y));
    else return vec3(0.0);
}

// Converts depth in [0,1] back to linear depth; adjust if you use non-linear depth
// Assuming linear depth was stored scaled to [0,1] range; else use your method
float LinearizeDepth(float depth) {
    const float near = 0.1;
    const float far = 100.0;
    return depth * (far - near) + near;
}

// Compute world position given probe position, direction, and depth
vec3 ReconstructPosition(vec3 dir, float linearDepth) {
    return probePosition.probePos + dir * linearDepth;
}

// Depth test - returns true if direction is unoccluded
bool IsVisible(vec3 dir, float surfaceDepth) {
    // Sample stored depth from depth cube map for this direction
    float storedDepth = texture(depthMap, vec4(dir,0)).r;
    if (storedDepth <= 0.0 || storedDepth >= 1.0) {
        // No valid geometry in this direction - treat as visible
        return true;
    }

    // Linearize depths if needed
    float sampledDepth = LinearizeDepth(storedDepth);
    float currentDepth = surfaceDepth;

    // Apply bias epsilon to avoid self-occlusion; tune as necessary
    const float depthBias = 0.001;

    // If current depth is significantly greater than stored, it means occluded
    if (currentDepth > sampledDepth + depthBias) {
        return false; // occluded
    }
    return true; // visible
}

void main() {
    // Reset SH coefficients
    for (int i = 0; i < NUM_COEFFS; ++i) {
        shBuffer.shCoeffs[i] = vec3(0.0);
    }

    // We'll sample directions using some stratified pattern on the sphere
    // For simplicity, sample uniformly over the sphere (e.g. cube map texels or uniform samples)
    // Here, we iterate over a sample count along a deterministic pattern

    // Accumulator for weight sums (for normalization)
    float weightSum = 0.0;

    // We'll sample uniformly distributed directions over the unit sphere
    // You may use better distributions (e.g., Hammersley) for improved results

    for (uint i = 0; i < NUM_SAMPLES; ++i) {
        // Generate direction sample
        // Simple: pick uniformly on sphere - use spherical coordinates
        float u = float(i) / float(NUM_SAMPLES);
        float v = fract(sin(float(i)*12.9898 + 78.233) * 43758.5453); // pseudo random in [0,1]

        // Spherical coords (theta from 0 to pi, phi from 0 to 2pi)
        float theta = acos(1.0 - 2.0 * u); // polar angle
        float phi = 2.0 * PI * v;          // azimuthal angle

        vec3 dir = vec3(
            sin(theta) * cos(phi),
            cos(theta),
            sin(theta) * sin(phi)
        );

        // Sample irradiance from irradiance cube map in direction
        vec3 radiance = texture(irradianceMap, dir).rgb;

        // Sample depth in that direction
        float sampledDepth = length(texture(depthMap, vec4(dir,0)).xyz);
        if (sampledDepth <= 0.1 || sampledDepth >= 100.0) {
            // No valid depth, assume visible
            // No occlusion
        } else {
            // Linearize
            float linearDepth = LinearizeDepth(sampledDepth);

            // Compute current depth as distance from probe in this dir to scene surface
            float currentDepth = length(dir) * linearDepth; // length(dir) = 1.0, so currentDepth = linearDepth

            // Depth test
            if (!IsVisible(dir, currentDepth)) {
                // Sample direction is occluded, skip contribution
                continue;
            }

            // Reconstruct sampled position
            //vec3 sampledPos = probePosition.probePos + dir * linearDepth;

            // Sample normal at this direction
            //vec3 sampledNormal = normalize(texture(normalMap, dir).xyz * 2.0 - 1.0);

            // Vector from sampled surface to probe
            //vec3 toProbe = normalize(probePosition.probePos - sampledPos);

            //const float normalFacingThreshold = 0.2; // tune this

            // Cull surfaces that face away from probe to avoid light bleeding
            //if (dot(sampledNormal, toProbe) < normalFacingThreshold) {
            //    continue;
            //}
        }

        // Apply cosine lobe weight (Lambert's law)
        float cosTheta = max(dir.y, 0.0); // IMPORTANT: Must be cosine with *local* normal = (0,1,0) for irradiance?
        // Actually, for environment lighting the normal is arbitrary; we can treat weight as 1/NUM_SAMPLES 
        // or weight by cosine if you want to simulate irradiance from hemisphere. Since you are capturing from all directions,
        // cosine weighting is needed when integrating irradiance over hemisphere aligned with surface normal.
        // But since the probe irradiance is over all directions, no cosine weight here; instead you only accumulate spherical harmonics weighted by 1/(4pi).

        // Differential solid angle for uniform samples: 4*PI / NUM_SAMPLES
        float dOmega = 4.0 * PI / float(NUM_SAMPLES);

        // Accumulate SH coefficients
        for (int c = 0; c < NUM_COEFFS; ++c) {
            vec3 basisVal = SHBasis(c, dir);
            // Radiance multiplied by basis * differential solid angle
            shBuffer.shCoeffs[c] += radiance * basisVal * dOmega;
        }

        weightSum += 1.0;
    }

    // Normalize SH coefficients by weightSum if needed (should be no more than NUM_SAMPLES)
    // but since we multiply by differential solid angle they are already weighted for integral over sphere.
    // Also you may want to apply a scaling factor to match your lighting shader's SH conventions

    // Done
}

Lighting Fragment Shader Excerpt:

#version 450
#extension GL_EXT_scalar_block_layout : enable

layout(scalar, set = 0, binding = 18) buffer SHBuffer {
    vec3 shCoeffs[];
} shBuffer;
layout(set = 0, binding = 19) uniform samplerCubeArray SHposMap;

layout(set = 0, binding = 20) uniform ProbePosition {
    vec3 probePos;

// shCoeffs is an array of vec3, length 9
// normal is the world space normal vector (unit length)
vec3 evalDiffuseIrradiance(vec3 normal) {
    normal = normalize(normal);
    float x = normal.x;
    float y = normal.y;
    float z = normal.z;

    float shBasis[9];
    shBasis[0] = 0.282095;         // 0: Y0^0
    shBasis[1] = 0.488603 * y;     // 1: Y1^-1
    shBasis[2] = 0.488603 * z;     // 2: Y1^0
    shBasis[3] = 0.488603 * x;     // 3: Y1^1
    shBasis[4] = 1.092548 * x * y; // 4: Y2^-2
    shBasis[5] = 1.092548 * y * z; // 5: Y2^-1
    shBasis[6] = 0.315392 * (3.0 * z * z - 1.0); // 6: Y2^0
    shBasis[7] = 1.092548 * x * z; // 7: Y2^1
    shBasis[8] = 0.546274 * (x * x - y * y); // 8: Y2^2

    vec3 irradiance = vec3(0.0);
    for (int i = 0; i < 9; ++i) {
        irradiance += shBuffer.shCoeffs[i] * shBasis[i];
    }

    // Clamp to zero or small epsilon to avoid negative lighting artifacts
    irradiance = max(irradiance, vec3(0.0));

    return irradiance;
}

float ComputeProbeOcclusion2(vec3 fragPosWorld, vec3 fragNormal)
{
    vec3 dir = normalize(fragPosWorld - probePosition.probePos);
    dir.y = -dir.y; // keep this if you need it for coordinate system alignment

    float fragDistance = length(fragPosWorld - probePosition.probePos);

    // Sample depth from cubemap [position stored in world space]
    float sampledDepth = length(texture(SHposMap, vec4(dir,0)).xyz - probePosition.probePos);

    const float depthBias = 0.002;
    float epsilon = 0.005;

    // Sample normal from cubemap at direction
    //vec3 sampledNormal = normalize(texture(SHnormalMap, dir).xyz);// * 2.0 - 1.0); // assuming normals stored in [0,1] range and remapped to [-1,1]

    // Vector from sampled position to the probe (points *towards* the probe)
    vec3 toProbe = normalize(probePosition.probePos - texture(SHposMap, vec4(dir,0)).xyz);
    toProbe.z = -toProbe.z;

    // Dot product test: is the surface facing the probe?
    float NdotV = dot(-toProbe, fragNormal);

    // Threshold to consider surface facing the probe—adjust tolerance as needed
    const float normalFacingThreshold = 0.001;

    // Compute occlusion only if surface faces probe
    float occlusion = 1.0;

    if (NdotV > normalFacingThreshold)
    {
        // If fragment farther than stored depth + bias => occluded
        occlusion = smoothstep(sampledDepth + depthBias + epsilon, sampledDepth + depthBias, fragDistance);
        occlusion = max(occlusion, 0.1);
    }
    else
    {
        // Surface is backfacing probe in this direction, so no occlusion contribution
        occlusion = 0.0;
    }

    return occlusion;
}

    vec3 irradiance = evalDiffuseIrradiance(fragNormal);
    float occlusion = ComputeProbeOcclusion2(fragPosition, fragNormal); // fragWorldPos must be passed from g-buffer
vec3 ambientDiffuse = irradiance * fragAlbedo * (1.0 - metallic) * kPi * occlusion;
} probePosition;