Skip to content

Commit 248ce25

Browse files
pezcodepezcode
committed
Add multiple scattering correction for dielectrics
1 parent 0f8e237 commit 248ce25

3 files changed

Lines changed: 85 additions & 80 deletions

File tree

‎src/Renderer/PBRShader.cpp‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -114,11 +114,10 @@ void PBRShader::generateAlbedoLUT()
114114
bgfx::blit(1, readbackTexture, 0, 0, albedoLUTTexture);
115115
116116
char* mem = new char[w * h * 4];
117-
bgfx::readTexture(readbackTexture, mem, 0);
117+
uint64_t frame = bgfx::readTexture(readbackTexture, mem, 0);
118118
119-
// wait two frames for result to be available and write it to a file
120-
bgfx::frame();
121-
bgfx::frame();
119+
// wait for result to be available and write it to a file
120+
while(bgfx::frame() < frame) { }
122121
123122
bx::FileWriter writer;
124123
if(writer.open(file, false, &err))

‎src/Renderer/Shaders/cs_multiple_scattering_lut.sc‎

Lines changed: 70 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,8 @@
66

77
// compute shader to calculate a lookup table for multiple scattering correction
88

9+
// TODO fix black pixels at grazing angles (NaN?)
10+
911
// Turquin. 2018. Practical multiple scattering compensation for microfacet models.
1012
// https://blog.selfshadow.com/publications/turquin/ms_comp_final.pdf
1113

@@ -16,8 +18,8 @@
1618
#define LUT_SIZE 32
1719
#define THREADS LUT_SIZE
1820

19-
// number of samples for approximating the albedo integral
20-
#define NUM_SAMPLES 2048
21+
// number of samples for approximating the integral
22+
#define NUM_SAMPLES 1024
2123

2224
// http://holger.dammertz.org/stuff/notes_HammersleyOnHemisphere.html
2325
vec2 hammersley_2d(uint i, uint N)
@@ -27,40 +29,38 @@ vec2 hammersley_2d(uint i, uint N)
2729
bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
2830
bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
2931
bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
30-
float ri = float(bits) * 2.3283064365386963e-10; // / uintBitsToFloat(0x100000000);
32+
float ri = float(bits) * 2.3283064365386963e-10; // / (1.0 << 32)
3133
return vec2(float(i) / float(N), ri);
3234
}
3335

34-
// map two random variables in [0;1] to a point on a hemisphere centered around N = y
36+
// map from [0;1] to a cosine weighted distribution on a hemisphere centered around N = y
3537
// http://holger.dammertz.org/stuff/notes_HammersleyOnHemisphere.html
36-
vec3 sampleHemisphere(vec2 Xi)
38+
vec3 sampleHemisphereCosine(vec2 Xi)
3739
{
3840
// turn Xi into spherical coordinates
3941
// azimuthal angle (360°)
4042
float phi = Xi.y * 2.0 * PI;
4143
// polar angle (90°)
4244
// 1-u to map the first sample (0) in the Hammersley sequence to 1=cos(0°)=up
43-
float cosTheta = 1.0 - Xi.x;
45+
float cosTheta = sqrt(1.0 - Xi.x);
4446
float sinTheta = sqrt(1.0 - cosTheta * cosTheta);
47+
4548
// to cartesian coordinates, y is up
4649
return vec3(
4750
cos(phi) * sinTheta,
4851
cosTheta,
4952
sin(phi) * sinTheta);
5053
}
5154

52-
// same thing, but importance sampled for the GGX NDF
53-
// return values is the half-way vector H
54-
vec3 importanceSampleGGX(vec2 Xi, float roughness)
55+
// importance sample the hemisphere for the GGX NDF instead
56+
// return value is the half-way vector H
57+
// a is perceptual roughness squared
58+
vec3 sampleGGX(vec2 Xi, float a)
5559
{
56-
float a = roughness * roughness;
57-
58-
// Sample in spherical coordinates
5960
float phi = 2.0 * PI * Xi.y;
6061
float cosTheta = sqrt((1.0 - Xi.x) / (1.0 + (a * a - 1.0) * Xi.x));
6162
float sinTheta = sqrt(1.0 - cosTheta * cosTheta);
6263

63-
// to cartesian coordinates, y is up
6464
return vec3(
6565
cos(phi) * sinTheta,
6666
cosTheta,
@@ -69,84 +69,74 @@ vec3 importanceSampleGGX(vec2 Xi, float roughness)
6969

7070
#define IMPORTANCE_SAMPLE_BRDF 1
7171

72-
// calculate the albedo for a perfectly reflective surface (F0 = 1)
73-
// = the irradiance reflected towards the eye position from uniform
74-
// lighting over the hemisphere
75-
float albedo(float NoV, float roughness)
72+
// calculate the directional albedo = the irradiance reflected towards the eye position
73+
// from uniform lighting over the hemisphere
74+
float albedo_specular(vec3 V, float NoV, PBRMaterial mat)
7675
{
77-
vec3 V;
78-
V.x = sqrt(1.0 - NoV * NoV); // sin
79-
V.y = NoV; // cos
80-
V.z = 0.0;
81-
8276
// N points straight upwards (y) for this integration
8377
const vec3 N = vec3(0.0, 1.0, 0.0);
8478

8579
float E = 0.0;
8680

87-
// fixed F0 -> perfectly reflective surface
88-
// only valid for metals, for dielectrics this needs to be a 3D LUT with F0 as a parameter
89-
// we could possibly fix F0 at 0.04 like GLTF does, and write a second channel in the LUT
90-
// how does this work with lerping with the metallic factor?
91-
const vec3 F0 = vec3_splat(1.0);
81+
// Monte-Carlo sampling for numerically integrating the directional albedo
82+
// E(V) = integral(brdf(V,L) * dot(N,L))
9283

93-
// Monte-Carlo sampling for numerically integrating the albedo over the hemisphere
94-
// a(V) = integral(brdf(V, L) * dot(N,L))
9584
for(uint i = 0; i < NUM_SAMPLES; i++)
9685
{
9786
// quasirandom values in [0;1]
9887
// is there a better distribution we can use?
9988
// see http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/
10089
vec2 Xi = hammersley_2d(i, NUM_SAMPLES);
10190

102-
#if IMPORTANCE_SAMPLE_BRDF
103-
10491
// sample microfacet direction
10592
// this returns H because we sample the NDF, which is the distribution of microfacets, which reflect around H
106-
// apparently the PDF needs the 4*VoH denominator (related to the Jacobian)
107-
// clarify!
108-
vec3 H = importanceSampleGGX(Xi, roughness);
93+
// then the PDF needs the 4*VoH denominator since we evaluate L (related to the Jacobian)
94+
vec3 H = sampleGGX(Xi, mat.a);
10995

11096
// get the light direction
11197
vec3 L = 2.0 * dot(V, H) * H - V;
11298

113-
#else
114-
115-
vec3 L = sampleHemisphere(Xi);
116-
vec3 H = normalize(V + L);
117-
118-
#endif
119-
12099
float NoL = saturate(dot(N, L));
121100
float NoH = saturate(dot(N, H));
122101
float VoH = saturate(dot(V, H));
123102

124-
// specular BRDF
125-
126-
float a = roughness * roughness;
127-
128-
float F = F_Schlick(VoH, F0).x;
129-
float VF = V_SmithGGXCorrelated(NoV, NoL, a);
130-
float D = D_GGX(NoH, a);
103+
// evaluate BRDF
131104

132-
#if IMPORTANCE_SAMPLE_BRDF
105+
float F = F_Schlick(VoH, mat.F0).x;
106+
float VF = V_SmithGGXCorrelated(NoV, NoL, mat.a);
107+
float D = D_GGX(NoH, mat.a);
133108

134109
//float Fr = F * VF * D;
135-
//float pdf = D * NoH / (4.0 * VoH);
136-
//E += Fr * NoL / pdf;
110+
//float inv_pdf = (4.0 * VoH) / (D * NoH);
111+
//E += Fr * NoL * inv_pdf;
137112
// -> D cancels out
138113

139-
E += F * VF * NoL * (4.0 * VoH) / NoH;
114+
E += F * VF * NoL * 4.0 * VoH / NoH;
115+
}
140116

141-
#else
117+
return E / float(NUM_SAMPLES);
118+
}
142119

143-
// this kind of converges at massive sample rates (1 million +)
144-
// but for low angles and high roughness the output is off
145-
// investigate!
146-
float pdf = 0.5 * INV_PI;
147-
E += F * VF * D * NoL / pdf;
120+
float albedo_diffuse(vec3 V, float NoV, PBRMaterial mat)
121+
{
122+
float E = 0.0;
123+
124+
for(uint i = 0; i < NUM_SAMPLES; i++)
125+
{
126+
vec2 Xi = hammersley_2d(i, NUM_SAMPLES);
127+
128+
vec3 L = sampleHemisphereCosine(Xi);
129+
vec3 H = normalize(V + L);
130+
131+
float VoH = saturate(dot(V, H));
148132

149-
#endif
133+
float F = F_Schlick(VoH, mat.F0).x;
134+
135+
// Fr = (1 - F) * C * (1/pi) * NoL
136+
// float inv_pdf = pi/NoL
137+
// -> 1/pi and NoL cancel out
138+
139+
E += (1.0 - F) * mat.diffuseColor.x;
150140
}
151141

152142
return E / float(NUM_SAMPLES);
@@ -160,9 +150,27 @@ void main()
160150

161151
float NoV = values.x;
162152
float roughness = values.y;
163-
float result = albedo(NoV, roughness);
164153

165-
result = saturate(result);
154+
vec3 V;
155+
V.x = sqrt(1.0 - NoV * NoV); // sin
156+
V.y = NoV; // cos
157+
V.z = 0.0;
158+
159+
PBRMaterial mat;
160+
// D3D compiler insists we initialize everything
161+
mat.normal = vec3(0.0, 1.0, 0.0);
162+
mat.occlusion = 0.0;
163+
mat.emissive = vec3_splat(0.0);
164+
mat.albedo = vec4_splat(1.0);
165+
mat.roughness = roughness;
166+
167+
mat.metallic = 1.0; // F0 = albedo -> perfectly reflective surface
168+
mat = pbrInitMaterial(mat);
169+
float albedo_metal = albedo_specular(V, NoV, mat);
170+
171+
mat.metallic = 0.0; // F0 = 0.04
172+
mat = pbrInitMaterial(mat);
173+
float albedo_dielectric = albedo_specular(V, NoV, mat) + albedo_diffuse(V, NoV, mat);
166174

167-
imageStore(i_texAlbedoLUT, coords, vec4(result, result, result, 1.0));
175+
imageStore(i_texAlbedoLUT, coords, vec4(albedo_metal, albedo_dielectric, 0.0, 1.0));
168176
}

‎src/Renderer/Shaders/pbr.sh‎

Lines changed: 12 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -161,6 +161,8 @@ PBRMaterial pbrInitMaterial(PBRMaterial mat)
161161
mat.F0 = mix(dielectricSpecular, mat.albedo.rgb, mat.metallic);
162162
// perceptual roughness to roughness
163163
mat.a = mat.roughness * mat.roughness;
164+
// prevent division by 0
165+
mat.a = max(mat.a, 0.01);
164166
165167
return mat;
166168
}
@@ -185,11 +187,7 @@ float specularAntiAliasing(vec3 N, float a)
185187
vec3 dndv = dFdy(N);
186188
float variance = SIGMA2 * (dot(dndu, dndu) + dot(dndv, dndv));
187189
float kernelRoughness2 = min(2.0 * variance, KAPPA);
188-
float filteredRoughness2 = saturate(a + kernelRoughness2);
189-
a = filteredRoughness2;
190-
191-
// Frostbite clamps roughness to 0.045 (0.045^2 = 0.002025)
192-
return max(a, 0.002025);
190+
return saturate(a + kernelRoughness2);
193191
}
194192

195193
#endif
@@ -266,20 +264,20 @@ float Fd_Lambert()
266264
// https://blog.selfshadow.com/publications/turquin/ms_comp_final.pdf
267265
vec3 multipleScatteringFactor(PBRMaterial mat, float NoV)
268266
{
269-
// E is the albedo for single scattering, ie. the total reflectance for a viewing direction
270-
float E = texture2D(s_texAlbedoLUT, vec2(NoV, mat.a)).x;
271-
vec3 factor = vec3_splat(1.0) + mat.F0 * (1.0/E - 1.0);
272-
273-
// TODO implement for dielectrics
274-
// requires an extra dimension in the LUT texture
267+
// Turquin approximates the multiple scattering portion of the BRDF using a scaled down version of the single scattering BRDF
268+
// That scale factor is E: the directional albedo for single scattering, ie. the total reflectance for a viewing direction
269+
vec2 E = texture2D(s_texAlbedoLUT, vec2(NoV, mat.a)).xy;
275270

276-
// for metals, the albedo value is calculated with F = 1
271+
// for metals, the albedo value is calculated with F = 1 (perfect reflection)
277272
// fresnel determines whether light is reflected or absorbed
273+
vec3 factorMetallic = vec3_splat(1.0) + mat.F0 * (1.0 / E.x - 1.0);
278274

279-
// for dielectrics, fresnel determines the ratio between two different albedos
275+
// for dielectrics, fresnel determines the ratio between specular and diffuse energy
280276
// so the albedo depends on F as a variable
277+
// however, dielectrics in GLTF have a fixed F0 of 0.04 so we can do this with a second LUT
278+
vec3 factorDielectric = vec3_splat(1.0 / E.y);
281279

282-
return mix(vec3_splat(1.0), factor, mat.metallic);
280+
return mix(factorDielectric, factorMetallic, mat.metallic);
283281
}
284282

285283
#endif

0 commit comments

Comments
 (0)