Spiral sampling over a solid angle

spiral_sampling.hlsl

// The vector around which the sampling solid-angle is centered.
// Does not matter for a full sphere, but does matter for anything smaller.
float3 central_direction = ??;
// [0, 1], where 0 is a single vector along the `central_direction` and 1 is a full sphere.
float angle_coefficient = ??;

const int NUM_SAMPLES = 100;
const float GOLDEN_RATIO_MULTIPLIED_BY_TWO_PI = 10.1664073846;

// Rotation matrix from Z-axis to the central_direction vector
// Source: https://math.stackexchange.com/a/476311
float3x3 rot;
{
    // central_direction must be normalized
    central_direction = normalize(central_direction);
    
    float3 axis = float3(0, 0, 1);
    float3 v = cross(axis, central_direction);
    float c = dot(axis, central_direction);

    float3x3 vx = {
        0.0, -v.z, v.y,
        v.z, 0.0, -v.x,
        -v.y, v.x, 0.0
    };

    float3x3 identity = {
        1.0, 0.0, 0.0,
        0.0, 1.0, 0.0,
        0.0, 0.0, 1.0,
    };

    rot =
        identity
        + vx
        + 1.0 / (1.0 + c) * mul(vx, vx);
}

for (int i = 0; i < NUM_SAMPLES; i++) {
    float phi = acos(1.0 - 2.0 * i / (float)NUM_SAMPLES); // [0, pi)
    phi *= angle_coefficient;
    float theta = GOLDEN_RATIO_MULTIPLIED_BY_TWO_PI * i;
    // `raw_direction` is centered around the Z-axis
    float3 raw_direction = float3(
        cos(theta) * sin(phi),
        sin(theta) * sin(phi),
        cos(phi)
    );
    // `sample_direction` is centered around the central_direction
    float3 sample_direction = mul(rot, raw_direction);
    
    // Do the actual sampling
    // ...
}