Sampling Visible GGX Normals with Spherical Caps

SampleVndf_GGX.cpp

// Helper function: sample the visible hemisphere from a spherical cap
vec3 SampleVndf_Hemisphere(vec2 u, vec3 wi)
{
	// sample a spherical cap in (-wi.z, 1]
	float phi = 2.0f * M_PI * u.x;
	float z = fma((1.0f - u.y), (1.0f + wi.z), -wi.z);
	float sinTheta = sqrt(clamp(1.0f - z * z, 0.0f, 1.0f));
	float x = sinTheta * cos(phi);
	float y = sinTheta * sin(phi);
	vec3 c = vec3(x, y, z);
	// compute halfway direction;
	vec3 h = c + wi;
	// return without normalization as this is done later (see line 25)
	return h;
}

// Sample the GGX VNDF
vec3 SampleVndf_GGX(vec2 u, vec3 wi, vec2 alpha)
{
    // warp to the hemisphere configuration
    vec3 wiStd = normalize(vec3(wi.xy * alpha, wi.z));
    // sample the hemisphere
    vec3 wmStd = SampleVndf_Hemisphere(u, wiStd);
    // warp back to the ellipsoid configuration
    vec3 wm = normalize(vec3(wmStd.xy * alpha, wmStd.z));
    // return final normal
    return wm;
}

#if 0 // alternative version from the two functions above: single function version 
vec3 SampleVndf_GGX(vec2 u, vec3 wi, vec2 alpha)
{
    // warp to the hemisphere configuration
    vec3 wiStd = normalize(vec3(wi.xy * alpha, wi.z));
    // sample a spherical cap in (-wi.z, 1]
    float phi = (2.0f * u.x - 1.0f) * M_PI;
    float z = fma((1.0f - u.y), (1.0f + wiStd.z), -wiStd.z);
    float sinTheta = sqrt(clamp(1.0f - z * z, 0.0f, 1.0f));
    float x = sinTheta * cos(phi);
    float y = sinTheta * sin(phi);
    vec3 c = vec3(x, y, z);
    // compute halfway direction as standard normal
    vec3 wmStd = c + wiStd;
    // warp back to the ellipsoid configuration
    vec3 wm = normalize(vec3(wmStd.xy * alpha, wmStd.z));
    // return final normal
    return wm;
}
#endif

#if 0 // world-space isotropic-only version
vec3 SampleVndf_GGX(vec2 u, vec3 wi, float alpha, vec3 n)
{
    // decompose the vector in parallel and perpendicular components
    vec3 wi_z = n * dot(wi, n);
    vec3 wi_xy = wi - wi_z;
    // warp to the hemisphere configuration
    vec3 wiStd = normalize(wi_z - alpha * wi_xy);
    // sample a spherical cap in (-wiStd.z, 1]
    float wiStd_z = dot(wiStd, n);
    float phi = (2.0f * u.x - 1.0f) * M_PI;
    float z = (1.0f - u.y) * (1.0f + wiStd_z) - wiStd_z;
    float sinTheta = sqrt(clamp(1.0f - z * z, 0.0f, 1.0f));
    float x = sinTheta * cos(phi);
    float y = sinTheta * sin(phi);
    vec3 cStd = vec3(x, y, z);
    // reflect sample to align with normal
    vec3 up = vec3(0, 0, 1);
    vec3 wr = n + up;
    vec3 c = dot(wr, cStd) * wr / wr.z - cStd;
    // compute halfway direction as standard normal
    vec3 wmStd = c + wiStd;
    vec3 wmStd_z = n * dot(n, wmStd);
    vec3 wmStd_xy = wmStd_z - wmStd;
    // warp back to the ellipsoid configuration
    vec3 wm = normalize(wmStd_z + alpha * wmStd_xy);
    // return final normal
    return wm;
}
// benefits: 
// - no need for moving to tangent space
// - it avoids the need for an orthonormal basis
// - it's (slightly) faster than the general version
#endif

Hi,

float z = 1 - u.y - u.y * wi.z;

Yes this is right.

Can we use this expanded term instead of fma without introducing precision problem？

I would certainly think so because fma implementations are generally less precise than the expanded form.
But don't necessarily take my word for it, feel free to experiment on your side :)

what would the pdf code looks like for this ?

The world-space isotropic-only version sometimes generated NaNs for me, to prevent this I added two checks which work for me:

vec3 SampleVndf_GGX(vec2 u, vec3 wi, float alpha, vec3 n)
{
+   // Dirac function for alpha = 0
+   if (alpha == 0.0f) return n;
    // decompose the vector in parallel and perpendicular components
    vec3 wi_z = n * dot(wi, n);
    vec3 wi_xy = wi - wi_z;
    // warp to the hemisphere configuration
    vec3 wiStd = normalize(wi_z - alpha * wi_xy);
    // sample a spherical cap in (-wiStd.z, 1]
    float wiStd_z = dot(wiStd, n);
    float phi = (2.0f * u.x - 1.0f) * M_PI;
    float z = (1.0f - u.y) * (1.0f + wiStd_z) - wiStd_z;
    float sinTheta = sqrt(clamp(1.0f - z * z, 0.0f, 1.0f));
    float x = sinTheta * cos(phi);
    float y = sinTheta * sin(phi);
    vec3 cStd = vec3(x, y, z);
    // reflect sample to align with normal
    vec3 up = vec3(0, 0, 1);
    vec3 wr = n + up;
-   vec3 c = dot(wr, cStd) * wr / wr.z - cStd;
+   // prevent division by zero
+   float wrz_safe = max(wr.z, 1e-6);
+   vec3 c = dot(wr, cStd) * wr / wrz_safe - cStd;
    // compute halfway direction as standard normal
    vec3 wmStd = c + wiStd;
    vec3 wmStd_z = n * dot(n, wmStd);
    vec3 wmStd_xy = wmStd_z - wmStd;
    // warp back to the ellipsoid configuration
    vec3 wm = normalize(wmStd_z + alpha * wmStd_xy);
    // return final normal
    return wm;
}

Are these correct?