Continue work on glitter

5 files changed, 201 insertions, 83 deletions

diff --git a/brdf.cginc b/brdf.cginc
index 54340aa..de222e7 100755
--- a/brdf.cginc
+++ b/brdf.cginc
@@ -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);
+}
+
diff --git a/math.cginc b/math.cginc
index 6bca086..55e1e35 100755
--- a/math.cginc
+++ b/math.cginc
@@ -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