diff options
Diffstat (limited to 'glitter.cginc')
| -rw-r--r-- | glitter.cginc | 152 |
1 files changed, 152 insertions, 0 deletions
diff --git a/glitter.cginc b/glitter.cginc new file mode 100644 index 0000000..39b735b --- /dev/null +++ b/glitter.cginc @@ -0,0 +1,152 @@ +#ifndef __GLITTER_INC +#define __GLITTER_INC + +#include "math.cginc" + +/* +@article{KPT:2025:Glinty, + title = {Evaluating and Sampling Glinty NDFs in Constant Time}, + author = {Kemppinen, Pauli and Paulin, LoÏs and Thonat, Théo and Thiery, Jean-Marc and Lehtinen, Jaakko and Boubekeur, Tamy}, + year = {2025}, + journal = {ACM Transactions on Graphics (Proc. SIGGRAPH Asia 2025)}, + volume = {44}, + number = {6}, + articleno = {255}, +} +*/ +// Ported from: https://www.shadertoy.com/view/tcdGDl + +// Lambert azimuthal equal area projection +float2 lambert(float3 v) { + return v.xy / sqrt(1 + v.z); +} + +// v is a microfacet normal that has been squished according to alpha, a +// roughness parameter. +float3 ndf_to_disk_ggx(float3 v, float alpha) { + // Map `v` onto a hemisphere. + float3 hemi = float3(v.xy / alpha, v.z); + float denom = dot(hemi, hemi); + // Project onto circle with equal area projection, and remap from [-1, 1] + // to [0, 1]. + float2 v_disk = lambert(normalize(hemi)) * 0.5 + 0.5; + float jacobian_determinant = 1.0 / (alpha * alpha * denom * denom); + return float3(v_disk, jacobian_determinant); +} + +// Computes (M^T M)^-1 +float2x2 inv_quadratic(float2x2 M) { + float D = determinant(M); + float A = dot(M[0] / D, M[0] / D); + float B = -dot(M[0] / D, M[1] / D); + float C = dot(M[1] / D, M[1] / D); + return float2x2(C, B, B, A); +} + +float2x2 uv_ellipsoid(float2x2 uv_J) { + float2x2 Q = inv_quadratic(transpose(uv_J)); + float tr = 0.5 * (Q[0][0] + Q[1][1]); + float D = sqrt(max(0.0, tr * tr - determinant(Q))); + float l1 = tr - D; + float l2 = tr + D; + float2 v1 = float2(l1 - Q[1][1], Q[0][1]); + float2 v2 = float2(Q[1][0], l2 - Q[0][0]); + float2 n = 1.f/sqrt(float2(l1, l2)); + return float2x2(normalize(v1)*n.x, normalize(v2)*n.y); +} + +float QueryLod(float2x2 uv_J, float filter_size) { + float s0 = length(uv_J[0]), s1 = length(uv_J[1]); + return log2(max(s0, s1) * filter_size) + pow(2.0, filter_size); +} + +float normal(float2x2 cov, float2 x) { + return exp(-.5 * dot(x, mul(inverse(cov), x))) / (sqrt(determinant(cov)) * 2.0 * PI); +} + +float Rand2D(float2 x, float2 y, float l, uint i) { + uint2 ux = asuint(x); + uint2 uy = asuint(y); + uint ul = asuint(l); + return hash22_fast((ux>>16|ux<<16) ^ uy ^ ul ^ (i*0x124u)); +} + +float Rand1D(float2 x, float2 y, float l, uint i) { + return Rand2D(x, y, l, i).x; +} + +// Bürmann series, see https://en.wikipedia.org/wiki/Error_function +float erf(float x) { + float e = exp(-x*x); + return sign(x) * 2.0 * sqrt((1.0 - e) / PI) * + (sqrt(PI) * 0.5 + 31./200. * e - 341.0/8000.0 * e * e); +} + +float cdf(float x, float mu, float sigma) { + return 0.5 + 0.5 * erf((x-mu)/(sigma*sqrt(2.0))); +} + +float integrate_interval(float x, float size, float mu, float stdev, float lower_limit, float upper_limit) { + return cdf(min(x+size, upper_limit), mu, stdev) - cdf(max(x-size, lower_limit), mu, stdev); +} + +float integrate_box(float2 x, float2 size, float2 mu, float2x2 sigma, float2 lower_limit, float2 upper_limit) { + return + integrate_interval(x.x, size.x, mu.x, sqrt(sigma[0][0]), lower_limit.x, upper_limit.x) * + integrate_interval(x.y, size.y, mu.y, sqrt(sigma[1][1]), lower_limit.y, upper_limit.y); +} + +float compensation(float2 x_a, float2x2 sigma_a, float res_a) { + float containing = integrate_box(0.5, 0.5, x_a, sigma_a, 0.0, 1.0); + float explicitly_evaluated = integrate_box(round(x_a*res_a)/res_a, 1.0/res_a, x_a, sigma_a, 0, 1); + return containing - explicitly_evaluated; +} + +float ndf(float3 h, float alpha, float glint_alpha, float2 uv, float2x2 uv_J, float N, float filter_size) { + float res = sqrt(N); + float2 x_s = uv; + float3 x_a_and_d = ndf_to_disk_ggx(h, alpha); + float2 x_a = x_a_and_d.xy; + float d = x_a_and_d.z; + + float lambda = QueryLod(res * uv_J, filter_size); + + float D_filter = .0; + + for(float m = .0; m<2.; m += 1.) { + float l = floor(lambda) + m; + + float w_lambda = 1. - abs(lambda - l); + float res_s = res * pow(2., -l); + float res_a = pow(2., l); + + float2x2 uv_J2 = filter_size * uv_J; + float2x2 sigma_s = uv_J2 * transpose(uv_J2); + + float2x2 sigma_a = d * pow(glint_alpha, 2.) * float2x2(1., .0, .0, 1.); + + float2 base_i_a = clamp(round(x_a * res_a), 1., res_a-1.); + for(int j_a = 0; j_a < 4; ++j_a) { + float2 i_a = base_i_a + float2(int2(j_a, j_a/2)%2)-.5; + + float2 base_i_s = round(x_s * res_s); + for(int j_s = 0; j_s < 4; ++j_s) { + float2 i_s = base_i_s + float2(int2(j_s, j_s/2)%2)-.5; + + float2 g_s = (i_s + Rand2D(i_s, i_a, l, 1u) - .5) / res_s; + float2 g_a = (i_a + Rand2D(i_s, i_a, l, 2u) - .5) / res_a; + + float r = Rand1D(i_s, i_a, l, 4u); + float roulette = smoothstep(max(.0, r-.1), min(1.0, r+.1), w_lambda); + + D_filter += roulette * normal(sigma_a, x_a - g_a) * normal(sigma_s, x_s - g_s) / N; + } + } + D_filter += w_lambda * compensation(x_a, sigma_a, res_a); + } + + return D_filter * d / PI; +} + +#endif // __GLITTER_INC + |