#ifndef __SSFD_INC
#define __SSFD_INC

#include "globals.cginc"

#if defined(_SSFD)
float ssfd(float2 uv, float scale, float max_fwidth, float2 uv_offset, texture3D noise)
{
  //float uv_fw = fwidth(uv.x) + fwidth(uv.y);
  // Original paper uses SVD instead of fwidth.
  float2x2 M = float2x2(ddx(uv), ddy(uv));
  float2x2 MtM = mul(transpose(M), M);
  float trace = MtM[0][0] + MtM[1][1];
  float det = determinant(MtM);
  // Calculate eigenvalues using quadratic formula.
  float tmp = sqrt(trace * trace - 4 * det);
  float e1 = (trace + tmp) * 0.5;
  float e2 = (trace - tmp) * 0.5;
  float2 singular_values = sqrt(float2(e1, e2));
  // Logic from original paper: the smaller eigenvalue corresponds to the
  // largest amount of stretching, so we use it to determine when to
  // subdivide.
  float uv_fw = singular_values.y;
  uv_fw *= scale;

  uint width, height, depth;
  noise.GetDimensions(width, height, depth);
  float bayer_res = sqrt(depth);

  // Suppose max_fwidth is 1.
  // uv_fw is 16. That means UV is changing a lot per pixel. That means we want to shrink the scale of the UV.
  // Factor is 16.
  // log_2(factor) is 4.
  // Divide original by 16.
  float fw_factor = max(uv_fw / max_fwidth, 1e-6);
  // Fractal transitions need to happen in octaves to match the Bayer
  // self-similarity used by the reference implementation.
  float fractal_level = log2(fw_factor);
  float fractal_level_floor = floor(fractal_level);
  float fractal_remainder = fractal_level - fractal_level_floor;

  float uv_scale = exp2(-fractal_level_floor);
  uv *= uv_scale;
  uv += uv_offset * uv_scale;

  float n_layers = depth;
  float min_sub_layer = max(1.0, 0.25 * n_layers);
  float max_sub_layer = n_layers;

  // Only the top 3/4 of the layer stack is fractally self-similar.
  float sub_layer = lerp(min_sub_layer, max_sub_layer, 1 - fractal_remainder);
  float uvw = (sub_layer - 0.5) / n_layers;

  float3 uv_3d = float3(uv, uvw);

  float dither = noise.SampleLevel(linear_repeat_s, uv_3d, 0);

  return dither;
}
#endif  // _SSFD

void apply_ssfd(v2f i, float2 uv, LightData l, float3 normal, inout float3 albedo) {
#if defined(_SSFD)
  float ssfd_mask = ssfd(uv, _SSFD_Scale, _SSFD_Max_Fwidth, 0, _SSFD_Noise);
  float ssfd_mask_fw = fwidth(ssfd_mask);
  // TODO I think anti aliasing is probably broken
  //float light_amount = 1.0 - (dot(l.indirect.diffuse_dominant_dir, normal) * 0.5 + 0.5);
  float light_amount = 0;
  float ssfd_threshold = saturate(light_amount + _SSFD_Threshold);
  ssfd_mask = smoothstep(ssfd_threshold - ssfd_mask_fw * 0.5, ssfd_threshold + ssfd_mask_fw * 0.5, ssfd_mask);
  albedo = lerp(albedo, _SSFD_Tint.rgb, ssfd_mask);
#endif
}

#endif  // __SSFD_INC