SPHQuad

struct SphQuad {
    vec2 lo,hi;
    vec2 loSq,hiSq;
    float z, zSq; //
    
    float b0, b1, b0sq, k; // misc precomputed constants
    float S; // solid angle of ’Q’
};
    
void SphQuadSample(in vec3 o, float u, float v, out vec3 w, out float pdf) {
    SphQuad squad;
    // compute rectangle coords in local reference system
    squad.lo = vec2(-1.0)-o.xy;
    squad.hi = vec2(1.0)-o.xy;
    squad.loSq = squad.lo*squad.lo;
    squad.hiSq = squad.hi*squad.hi;
    squad.z = -abs(o.z);
    squad.zSq = o.z*o.z;
    // compute normals to edges
    float n0 = squad.lo.y*inversesqrt(squad.loSq.y+squad.zSq);
    float n1 = squad.hi.x*inversesqrt(squad.hiSq.x+squad.zSq);
    float n2 = squad.hi.y*inversesqrt(squad.hiSq.y+squad.zSq);
    float n3 = squad.lo.x*inversesqrt(squad.loSq.x+squad.zSq);
    // compute internal angles (gamma_i)
    float g0 = acos(n0*n1); // A/sqrt(A*A+Z*Z)*B/sqrt(B*B+Z*Z)
    float g1 = acos(-n1*n2);
    float g2 = acos(n2*n3);
    float g3 = acos(-n3*n0);
    // compute predefined constants
    squad.b0 = -n0;
    squad.b1 = n2;
    squad.b0sq = squad.loSq.y/(squad.loSq.y+squad.zSq);
    squad.k = 2.0*PI - g2 - g3;
    // compute solid angle from internal angles
    squad.S = g0 + g1 - squad.k;

    
    // 1. compute ’cu’
    float au = u * squad.S + squad.k;
    float fu = (cos(au) * squad.b0 - squad.b1) / sin(au);
    float cu = inversesqrt(fu*fu + squad.b0sq) * (fu>0. ? +1. : -1.);
    ///cu = clamp(cu, -1., 1.); // avoid NaNs
    // 2. compute ’xu’
    float xu = -(cu * squad.z)*inversesqrt(1. - cu*cu);
    ///xu = clamp(xu, squad.lo.x, squad.hi.x); // avoid Infs
    // 3. compute ’yv’
    float d2 = xu*xu + squad.zSq;
    float h0 = squad.lo.y*inversesqrt(d2 + squad.loSq.y);
    float h1 = squad.hi.y*inversesqrt(d2 + squad.hiSq.y);
    float hv = h0 + v * (h1-h0), hv2 = hv*hv;
    float yv = (hv2 < 1.-EPSILON) ? hv*inversesqrt((1.-hv2)/d2) : squad.hi.y;

    // TODO: transform (xu,yv,z) to world coords and compensate for solid angle
    w = normalize(vec3(xu,yv,squad.z));
    pdf = 1. / squad.S;
}