diff options
| -rwxr-xr-x | brdf.cginc | 5 | ||||
| -rw-r--r-- | glitter.cginc | 152 | ||||
| -rw-r--r-- | glitter.slang | 74 | ||||
| -rw-r--r-- | hashwithoutsine.cginc | 35 | ||||
| -rwxr-xr-x | math.cginc | 18 |
5 files changed, 201 insertions, 83 deletions
@@ -1,11 +1,12 @@ #ifndef __BRDF_INC #define __BRDF_INC -#include "pema99.cginc" -#include "pbr.cginc" +#include "glitter.cginc" #include "lighting.cginc" #include "lysenko.cginc" #include "math.cginc" +#include "pema99.cginc" +#include "pbr.cginc" float pow5(float x) { float x2 = x * x; 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 + diff --git a/glitter.slang b/glitter.slang deleted file mode 100644 index 5442a0f..0000000 --- a/glitter.slang +++ /dev/null @@ -1,74 +0,0 @@ -#ifndef __GLITTER_INC -#define __GLITTER_INC - -#define PI 3.1415926535897932384626433832795028841971 - -float gaussian(float3 x, float3 mu, float sigma) { - float3 x_minus_mu = x - mu; - return exp(-dot(x_minus_mu, x_minus_mu) / (2 * sigma * sigma)) / (sigma * sqrt(2 * PI)); -} - -float det(float2x2 m) { - return m[0][0] * m[1][1] - m[0][1] * m[1][0]; -} - -float dirac(float3 x) { - // Ternary operator short-circuits in slang >:( - return any(abs(x) < 1e-4) ? 1e4 : 0; -} - -// Maps a point on the unit hemisphere to an inscribed circle. -float2 hemisphere_to_inscribed_circle(float3 p) { - // TODO - return 0; -} - -float3 rand3_hash(float3 p) -{ - // Improved Murmurhash3 by Squirrel Eiserloh (GDC 2017) - p = float3(dot(p, float3(127.1, 311.7, 74.7)), - dot(p, float3(269.5, 183.3, 246.1)), - dot(p, float3(113.5, 271.9, 124.6))); - return -1.0 + 2.0 * frac(sin(p) * 43758.5453123); -} - -public float D_Kemppinen( - float3 worldPos, - float3 H_tan, // halfway vector in tangent space - float amount, - float roughness) -{ - float2x2 J_T__omega_h = 0; // TODO - float det_J_T__omega_h = det(J_T__omega_h); - - float sum = 0; - // x is the current world space position. - float3 x = worldPos; - // omega_h is the halfway vector *in TBN coordinates.* - float3 omega_h = H_tan; - float3 T_omega_h = float3( - hemisphere_to_inscribed_circle(omega_h), - 0); - uint K_p = 1; // TODO - for (uint i = 0; i < K_p; ++i) { - // x_i is the position of the facet. - float3 x_i = worldPos; // TODO randomize position - // mu_i is the orientation of the facet. - float3 mu_i = abs(rand3_hash(x_i)); - - float a_star = roughness; - - float cur_gauss = gaussian( - T_omega_h, - mu_i, - a_star * sqrt(det_J_T__omega_h)); - - sum += cur_gauss * dirac(x_i - x); - } - - float scale = det_J_T__omega_h / K_p; - return sum * scale; -} - -#endif // __GLITTER_INC - diff --git a/hashwithoutsine.cginc b/hashwithoutsine.cginc new file mode 100644 index 0000000..78c6489 --- /dev/null +++ b/hashwithoutsine.cginc @@ -0,0 +1,35 @@ +/* Copyright (c)2014 David Hoskins. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE.*/ + +// Source: +// https://www.shadertoy.com/view/4djSRW + +float2 hash22_fast(float2 p) { + float3 p3 = frac(p.xyx * float3(0.1031, 0.1030, 0.0973)); + p3 += dot(p3, p3.yzx + 33.33); + return frac((p3.xx+p3.yz)*p3.zy); +} + +float3 hash33_fast(float3 p) { + p = frac(p * float3(0.1031, 0.1030, 0.0973)); + p += dot(p, p.yxz + 33.33); + return frac((p.xxy + p.yxx) * p.zyx); +} + @@ -1,6 +1,8 @@ #ifndef __MATH_INC #define __MATH_INC +#include "hashwithoutsine.cginc" + #define PI 3.14159265358979323846264f #define TAU (2.0f * PI) #define HALF_PI (PI * 0.5f) @@ -115,13 +117,6 @@ float4 alpha_blend(float4 front, float4 back) { return float4(front.rgb * front.a + back.rgb * (1 - front.a), front.a + back.a * (1 - front.a)); } -// 3 in 3 out hash. Based on the "hashwithoutsine" family. -float3 hash33_fast(float3 p) { - p = frac(p * float3(0.1031, 0.1030, 0.0973)); - p += dot(p, p.yxz + 33.33); - return frac((p.xxy + p.yxx) * p.zyx); -} - // O'Neill-style PCG 32-bit output permutation (RXS-M-XS). uint pcg32(uint input) { uint state = input * 747796405u + 2891336453u; @@ -263,4 +258,13 @@ float3 srgb_to_linear(float3 srgb_color) { return lerp(lo, hi, step(0.04045f, srgb_color)); } +float2x2 inverse(float2x2 m) { + float det = (m[0][0] * m[1][1]) - (m[0][1] * m[1][0]); + + return float2x2( + m[1][1], -m[0][1], + -m[1][0], m[0][0] + ) / det; +} + #endif // __MATH_INC |