As promised here's some code. I'm still not sure about some details (in particular those related to importance sampling) and it'd be great if other people could spot errors and add more variants to the code, like different distribution functions.
void GenerateFGLookupTexture(Graphics::Texture* texture, RenderingDevice& device)
{
Graphics::TextureDesc textureDesc;
device.GetDesc(texture, textureDesc);
assert(textureDesc.format == Graphics::Format::G16R16);
assert(textureDesc.levels == 1);
const float deltaNdotV = 1.f / static_cast<float>(textureDesc.width);
const float deltaRoughness = 1.f / static_cast<float>(textureDesc.height);
// Lock texture
const Graphics::MappedData mappedData = device.Map(texture, 0, Graphics::MapType::Write, 0);
uint8* const RESTRICT dataPtr = reinterpret_cast<uint8*>(mappedData.data);
const uint numSamples = 512;
float roughness = 0.f;
for (uint v = 0; v < textureDesc.height; v++)
{
float NdotV = 0.f;
uint32* dst = reinterpret_cast<uint32*>(dataPtr + v * mappedData.rowPitch);
for (uint u = 0; u < textureDesc.width; u++, ++dst)
{
float a = 0;
float b = 0;
CalculateFGCoeff(NdotV, roughness, numSamples, a, b);
assert(a <= 1.f);
assert(b <= 1.f);
*dst = (PACKINTOSHORT_0TO1(b) << 16) | PACKINTOSHORT_0TO1(a);
NdotV += deltaNdotV;
}
roughness += deltaRoughness;
}
device.Unmap(texture, 0);
}
void PlaneHammersley(float& x, float& y, int k, int n)
{
float u = 0;
float p = 0.5f;
// FIXME Optimize by removing conditional
for (int kk = k; kk; p *= 0.5f, kk /= 2)
{
if (kk & 1)
{
u += p;
}
}
x = u;
y = (k + 0.5f) / n;
}
vec3 ImportanceSampleGGX(float x, float y, float a4)
{
const float PI = 3.1415926535897932384626433832795028841971693993751f;
// Convert uniform random variables x, y to a sample direction
const float phi = 2 * PI * x;
const float cosTheta = std::sqrt( (1 - y) / ( 1 + (a4 - 1) * y) );
const float sinTheta = std::sqrt(1 - cosTheta * cosTheta);
// Convert direction to cartesian coordinates
const vec3 H(sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta);
//D = a2 / (PI * std::pow(cosTheta * (a2 - 1) + 1, 2));
return H;
}
// Reference: GPU Gems 3 - GPU Based Importance Sampling
//s2012_pbs_physics_math_slides.pdf
vec3 ImportanceSampleBlinn(float x, float y, float specularPower)
{
const float PI = 3.1415926535897932384626433832795028841971693993751f;
// Convert uniform random variables x, y to a sample direction
const float phi = 2 * PI * x;
const float cosTheta = std::pow(y, 1.f / (specularPower + 1));
const float sinTheta = std::sqrt(1 - cosTheta * cosTheta);
// Convert direction to cartesian coordinates
const vec3 H(sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta);
//D = (specularPower + 2) / (2 * PI) * std::pow(cosTheta, specularPower);
return H;
}
float G_Schlick(float k, float NdotV, float NdotL)
{
return (NdotV * NdotL) / ( (NdotV * (1 - k) + k) * (NdotL * (1 - k) + k) );
}
void CalculateFGCoeff(float NoV, float roughness, uint numSamples, float& a, float& b)
{
// Work in a coordinate system where normal = vec3(0,0,1)
#define GGX 0
#define BLINN 1
#define G_SCHLICK 0
// Build view vector
const vec3 V(std::sqrt(1.0f - NoV * NoV), 0, NoV);
#if BLINN
const float blinnSpecularPower = std::pow(2.f, 13.f * roughness);
#elif GGX
const float GXX_a4 = std::pow(roughness, 4);
#endif
#if G_SCHLICK
const float G_k = std::pow(roughness + 1, 2.f) / 8.f;
#endif
a = 0;
b = 0;
for (uint i = 0; i < numSamples; ++i)
{
float x, y;
PlaneHammersley(x, y, i, numSamples);
// Microfacet specular model:
// f = D*G*F / (4*NoL*NoV)
// V = G / (NoL*NoV)
// Importance-based sampling:
// f / pdf
// Calculate random half vector based on roughness
#if BLINN
const vec3 H = ImportanceSampleBlinn(x, y, blinnSpecularPower);
// D and pdfH cancel each other so just set to 1
const float D = 1;
const float pdfH = D;
#elif GXX
const vec3 H = ImportanceSampleGGX(x, y, GXX_a4);
// D and pdfH cancel each other so just set to 1
const float D = 1;
const float pdfH = D;
#endif
// Calculate light direction
const vec3 L = 2 * dot( V, H ) * H - V;
const float NoL = saturate( L.z );
const float VoH = saturate( dot( V, H ) );
const float NoH = saturate( H.z );
const float LoH = saturate( dot( L, H ) ); // = VoH
// Convert pdf(H) to pdf(L)
// Reference: Torrance and Sparrow
// http://graphics.stanford.edu/~boulos/papers/brdftog.pdf
const float pdfL = pdfH / (4 * LoH);
if (NoL > 0)
{
#if G_SCHLICK
// FIXME NoV cancel out
const float G = G_Schlick(G_k, NoV, NoL);
const float V = G / (NoL * NoV);
#elif
const float V = 1;
#endif
const float G_Vis = D * V / 4 / pdfL;
const float Fc = std::pow(1 - VoH, 5.f);
// FIXME : NoL ? Part of the lighting eq but gives dark reflections at grazing angles. Need a better BRDF probably
a += (1 - Fc) * G_Vis * NoL;
b += Fc * G_Vis * NoL;
}
}
a /= numSamples;
b /= numSamples;
}
}
