#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);
}

// Rebuild GLSL mat2 column semantics explicitly in HLSL.
float2 mat2_col(float2x2 m, uint i) {
  return float2(m[0][i], m[1][i]);
}

float2x2 mat2_from_cols(float2 c0, float2 c1) {
  return float2x2(c0.x, c1.x,
                  c0.y, c1.y);
}

// 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);
	float2 c0 = mat2_col(M, 0) / D;
	float2 c1 = mat2_col(M, 1) / D;
	float A = dot(c0, c0);
	float B = -dot(c0, c1);
	float C = dot(c1, c1);
	return mat2_from_cols(float2(C, B), float2(B, A));
}

float2x2 uv_ellipsoid(float2x2 uv_J) {
	float2x2 Q = inv_quadratic(transpose(uv_J));
	float2 q0 = mat2_col(Q, 0);
	float2 q1 = mat2_col(Q, 1);
	float tr = 0.5 * (q0.x + q1.y);
	float  D = sqrt(max(0.0, tr * tr - determinant(Q)));
	float l1 = tr - D;
	float l2 = tr + D;
	float2 v1 = float2(l1 - q1.y, q0.y);
	float2 v2 = float2(q1.x, l2 - q0.x);
	float2 n = 1.f/sqrt(float2(l1, l2));
	return mat2_from_cols(normalize(v1) * n.x, normalize(v2) * n.y);
}

float QueryLod(float2x2 uv_J, float filter_size) {
    float s0 = length(mat2_col(uv_J, 0));
    float s1 = length(mat2_col(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);
}

uint2 shuffle(uint2 v) {
    v = v * 1664525u + 1013904223u;
	v.x += v.y * 1664525u;
	v.y += v.x * 1664525u;

	v = v ^ (v>>16u);

	v.x += v.y * 1664525u;
	v.y += v.x * 1664525u;
	v = v ^ (v>>16u);
    return v;
}

float2 rand(uint2 v) {
	return float2(shuffle(v)) * pow(0.5, 32.0);
}

float2 Rand2D(float2 x, float2 y, float l, uint i) {
	uint2 ux = asuint(x);
	uint2 uy = asuint(y);
	uint  ul = asuint(l);
	return rand((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 D_Kemppinen(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 = mul(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