#ifndef __APERIODIC_TILING_INC
#define __APERIODIC_TILING_INC

// 5D cut-and-project for Penrose tiling.
// References:
//   https://www.shadertoy.com/view/XccXW8 (public domain)
//   https://gglouser.github.io/cut-and-project-tiling/docs/intro.html#how-it-works
//   https://www.youtube.com/watch?v=hwMAOFb6yvA

#include "globals.cginc"
#include "math.cginc"

#if defined(_APERIODIC_TILING)

static const float M5 = sqrt(2.0 / 5.0);
static const float APERIODIC_FILTER_THRESHOLD = 0.5;
static const float APERIODIC_FAR_THRESHOLD = 1.0;
static const float APERIODIC_FILTER_BLEND_WIDTH = 0.5;
static const float APERIODIC_SMOOTHSTEP_SQ_MEAN = 13.0 / 35.0;

static const float4 basis_u5_03 = M5 * float4(
  cos(0 * TAU / 10),
  cos(1 * TAU / 10),
  cos(2 * TAU / 10),
  cos(3 * TAU / 10));
static const float basis_u5_44 = M5 * cos(4 * TAU / 10);

static const float4 basis_v5_03 = M5 * float4(
  sin(0 * TAU / 10),
  sin(1 * TAU / 10),
  sin(2 * TAU / 10),
  sin(3 * TAU / 10));
static const float basis_v5_44 = M5 * sin(4 * TAU / 10);

static const float2 aperiodic_tile_offsets[4] = {
  float2(0.5, 0.5),
  float2(0.5, -0.5),
  float2(-0.5, 0.5),
  float2(-0.5, -0.5)
};

static const float4 aperiodic_face_a03[10] = {
  float4(1, 0, 0, 0),
  float4(1, 0, 0, 0),
  float4(1, 0, 0, 0),
  float4(1, 0, 0, 0),
  float4(0, 1, 0, 0),
  float4(0, 1, 0, 0),
  float4(0, 1, 0, 0),
  float4(0, 0, 1, 0),
  float4(0, 0, 1, 0),
  float4(0, 0, 0, 1)
};

static const float aperiodic_face_a44[10] = {
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0
};

static const float4 aperiodic_face_b03[10] = {
  float4(0, 1, 0, 0),
  float4(0, 0, 1, 0),
  float4(0, 0, 0, 1),
  float4(0, 0, 0, 0),
  float4(0, 0, 1, 0),
  float4(0, 0, 0, 1),
  float4(0, 0, 0, 0),
  float4(0, 0, 0, 1),
  float4(0, 0, 0, 0),
  float4(0, 0, 0, 0)
};

static const float aperiodic_face_b44[10] = {
  0, 0, 0, 1, 0, 0, 1, 0, 1, 1
};

static const float4 aperiodic_face_c03[10] = {
  float4(0, 0, 1, 1),
  float4(0, 1, 0, 1),
  float4(0, 1, 1, 0),
  float4(0, 1, 1, 1),
  float4(1, 0, 0, 1),
  float4(1, 0, 1, 0),
  float4(1, 0, 1, 1),
  float4(1, 1, 0, 0),
  float4(1, 1, 0, 1),
  float4(1, 1, 1, 0)
};

static const float aperiodic_face_c44[10] = {
  1, 1, 1, 0, 1, 1, 0, 1, 0, 0
};

// Area weights for the 10 face orientations. Adjacent axis pairs are the thin
// rhombs, separated pairs are the thick rhombs. These weights sum to 1.
static const float aperiodic_face_weights[10] = {
  0.0552786404500042,
  0.1447213595499958,
  0.1447213595499958,
  0.0552786404500042,
  0.0552786404500042,
  0.1447213595499958,
  0.1447213595499958,
  0.0552786404500042,
  0.1447213595499958,
  0.0552786404500042
};

// Precomputed per-face barycentric transforms. For each fixed face
// orientation, this is inverse(float2x2(proj5(a), proj5(b))).
static const float2x2 aperiodic_face_matrices[10] = {
  float2x2(1.5811388300841898, -2.1762508994828216, -0.0, 2.6899940478558295),
  float2x2(1.5811388300841895, -0.5137431483730078, -0.0, 1.6625077511098139),
  float2x2(1.5811388300841895, 0.5137431483730076, -0.0, 1.6625077511098136),
  float2x2(1.5811388300841895, 2.176250899482821, -0.0, 2.6899940478558286),
  float2x2(2.558336368008464, -0.8312538755549072, -1.58113883008419, 2.176250899482822),
  float2x2(1.5811388300841895, 0.5137431483730076, -0.9771975379242739, 1.3449970239279145),
  float2x2(0.977197537924274, 1.3449970239279145, -0.9771975379242739, 1.3449970239279145),
  float2x2(2.558336368008464, 0.8312538755549067, -2.5583363680084634, 0.8312538755549069),
  float2x2(0.9771975379242742, 1.3449970239279148, -1.5811388300841895, 0.5137431483730078),
  float2x2(1.5811388300841898, 2.176250899482821, -2.5583363680084634, -0.8312538755549066)
};

float dot5(float4 a03, float a44, float4 b03, float b44) {
  return dot(a03, b03) + a44 * b44;
}

float2 proj5(float4 p03, float p44) {
  return float2(
    dot5(p03, p44, basis_u5_03, basis_u5_44),
    dot5(p03, p44, basis_v5_03, basis_v5_44));
}

float3 aperiodic_face_color(int face_id) {
  if (face_id == 0) return _Aperiodic_Tiling_Color_0.rgb;
  if (face_id == 1) return _Aperiodic_Tiling_Color_1.rgb;
  if (face_id == 2) return _Aperiodic_Tiling_Color_2.rgb;
  if (face_id == 3) return _Aperiodic_Tiling_Color_3.rgb;
  if (face_id == 4) return _Aperiodic_Tiling_Color_4.rgb;
  if (face_id == 5) return _Aperiodic_Tiling_Color_5.rgb;
  if (face_id == 6) return _Aperiodic_Tiling_Color_6.rgb;
  if (face_id == 7) return _Aperiodic_Tiling_Color_7.rgb;
  if (face_id == 8) return _Aperiodic_Tiling_Color_8.rgb;
  return _Aperiodic_Tiling_Color_9.rgb;
}

float2 aperiodic_tile_barycentric(float2 uv_m, float2x2 m, float2 s, float2 q,
                                  float4 c03, float c44) {
  float4 shift03 = c03 * round(q.x * basis_u5_03 + q.y * basis_v5_03);
  float shift44 = c44 * round(q.x * basis_u5_44 + q.y * basis_v5_44);
  return uv_m - mul(m, proj5(shift03, shift44)) - s;
}

float interval_box_coverage(float center, float half_extent, float radius) {
  if (half_extent <= 1e-5) {
    return abs(center) <= radius ? 1.0 : 0.0;
  }

  float lo = center - half_extent;
  float hi = center + half_extent;
  float overlap = max(min(hi, radius) - max(lo, -radius), 0.0);
  return overlap / (2.0 * half_extent);
}

float square_box_coverage(float2 barycentric, float2 half_extents, float radius) {
  return interval_box_coverage(barycentric.x, half_extents.x, radius) *
         interval_box_coverage(barycentric.y, half_extents.y, radius);
}

#if defined(_APERIODIC_TILING_NORMALS)
float3 aperiodic_tiling_normal(float2 barycentric) {
  float bevel_width = min(_Aperiodic_Tiling_Normal_Thickness, 0.5);
  if (bevel_width <= 1e-5) {
    return float3(0.0, 0.0, 1.0);
  }

  float flat_limit = 0.5 - bevel_width;
  float2 edge_factor = smoothstep(flat_limit, 0.5, abs(barycentric));
  float2 xy = sign(barycentric) * edge_factor * _Aperiodic_Tiling_Normal_Strength;
  return normalize(float3(xy, 1.0));
}
#endif  // _APERIODIC_TILING_NORMALS

struct AperiodicPointSample {
  float distance_to_edge;
  float2 barycentric;
  int face_id;
};

struct AperiodicFootprintRange {
  float min_extent;
  float max_extent;
};

void aperiodic_accumulate_point_orientation(float2 uv, float4 p03, float p44, int face_id,
                                            inout AperiodicPointSample best) {
  float2x2 m = aperiodic_face_matrices[face_id];
  float4 a03 = aperiodic_face_a03[face_id];
  float a44 = aperiodic_face_a44[face_id];
  float4 b03 = aperiodic_face_b03[face_id];
  float b44 = aperiodic_face_b44[face_id];
  float4 c03 = aperiodic_face_c03[face_id];
  float c44 = aperiodic_face_c44[face_id];

  float2 r = round(float2(dot5(p03, p44, a03, a44), dot5(p03, p44, b03, b44)));
  float2 uv_m = mul(m, uv);

  [unroll]
  for (int candidate_id = 0; candidate_id < 4; ++candidate_id) {
    float2 s = r + aperiodic_tile_offsets[candidate_id];
    float2 q = mul(s, m);
    float2 barycentric = aperiodic_tile_barycentric(uv_m, m, s, q, c03, c44);
    float distance_to_edge = 0.5 - max(abs(barycentric.x), abs(barycentric.y));
    bool better = distance_to_edge > best.distance_to_edge;
    best.distance_to_edge = better ? distance_to_edge : best.distance_to_edge;
    best.barycentric = better ? barycentric : best.barycentric;
    best.face_id = better ? face_id : best.face_id;
  }
}

void sample_aperiodic_tiling_point(float2 uv, float4 p03, float p44, out float3 albedo,
                                   out float3 tiling_normal_tangent) {
  AperiodicPointSample best;
  best.distance_to_edge = -1e10;
  best.barycentric = 0.0;
  best.face_id = 0;

  [unroll]
  for (int face_id = 0; face_id < 10; ++face_id) {
    aperiodic_accumulate_point_orientation(uv, p03, p44, face_id, best);
  }

  float edge_width = min(_Aperiodic_Tiling_Edge_Thickness, 0.5);
  float edge_sd = best.distance_to_edge - edge_width;
  float edge_sd_aa = max(abs(fwidth(edge_sd)), 1e-4);
  float edge_mask = smoothstep(-edge_sd_aa * 0.5, edge_sd_aa * 0.5, edge_sd);
  albedo = lerp(_Aperiodic_Tiling_Edge_Color.rgb, aperiodic_face_color(best.face_id), edge_mask);

#if defined(_APERIODIC_TILING_NORMALS)
  float3 tile_normal = aperiodic_tiling_normal(best.barycentric);
  tiling_normal_tangent = normalize(float3(tile_normal.xy * edge_mask, 1.0));
#else
  tiling_normal_tangent = 0.0;
#endif
}

AperiodicFootprintRange aperiodic_footprint_range(float2 uv_ddx, float2 uv_ddy) {
  AperiodicFootprintRange range;
  range.min_extent = 1e10;
  range.max_extent = 0.0;

  [unroll]
  for (int face_id = 0; face_id < 10; ++face_id) {
    float2x2 m = aperiodic_face_matrices[face_id];
    float2 dbdx = mul(m, uv_ddx);
    float2 dbdy = mul(m, uv_ddy);
    float face_extent = max(abs(dbdx.x) + abs(dbdy.x), abs(dbdx.y) + abs(dbdy.y)) * 0.5;
    range.min_extent = min(range.min_extent, face_extent);
    range.max_extent = max(range.max_extent, face_extent);
  }

  return range;
}

void aperiodic_accumulate_filtered_orientation(float2 uv, float4 p03, float p44,
                                               float2 uv_ddx, float2 uv_ddy, int face_id,
                                               float inner_radius, inout float3 face_color_sum,
                                               inout float face_sum) {
  float2x2 m = aperiodic_face_matrices[face_id];
  float4 a03 = aperiodic_face_a03[face_id];
  float a44 = aperiodic_face_a44[face_id];
  float4 b03 = aperiodic_face_b03[face_id];
  float b44 = aperiodic_face_b44[face_id];
  float4 c03 = aperiodic_face_c03[face_id];
  float c44 = aperiodic_face_c44[face_id];

  float2 r = round(float2(dot5(p03, p44, a03, a44), dot5(p03, p44, b03, b44)));
  float2 uv_m = mul(m, uv);
  float2 dbdx = mul(m, uv_ddx);
  float2 dbdy = mul(m, uv_ddy);
  float2 half_extents = 0.5 * (abs(dbdx) + abs(dbdy));
  float2 inner_limit = half_extents + inner_radius;
  float3 face_color = aperiodic_face_color(face_id);

  [unroll]
  for (int candidate_id = 0; candidate_id < 4; ++candidate_id) {
    float2 s = r + aperiodic_tile_offsets[candidate_id];
    float2 q = mul(s, m);
    float2 barycentric = aperiodic_tile_barycentric(uv_m, m, s, q, c03, c44);

    if (abs(barycentric.x) >= inner_limit.x || abs(barycentric.y) >= inner_limit.y) {
      continue;
    }

    float inner = square_box_coverage(barycentric, half_extents, inner_radius);
    face_color_sum += face_color * inner;
    face_sum += inner;
  }
}

void sample_aperiodic_tiling_filtered(float2 uv, float4 p03, float p44, float2 uv_ddx,
                                      float2 uv_ddy, out float3 albedo,
                                      out float3 tiling_normal_tangent) {
  float inner_radius = max(0.5 - min(_Aperiodic_Tiling_Edge_Thickness, 0.5), 0.0);
  float3 face_color_sum = 0.0;
  float face_sum = 0.0;

  [unroll]
  for (int face_id = 0; face_id < 10; ++face_id) {
    aperiodic_accumulate_filtered_orientation(
        uv, p03, p44, uv_ddx, uv_ddy, face_id, inner_radius, face_color_sum, face_sum);
  }

  float edge_weight = saturate(1.0 - face_sum);
  albedo = face_color_sum + _Aperiodic_Tiling_Edge_Color.rgb * edge_weight;

#if defined(_APERIODIC_TILING_NORMALS)
  // In the filtered regime the bevel is subpixel detail, so keep the normal flat.
  tiling_normal_tangent = float3(0.0, 0.0, 1.0);
#else
  tiling_normal_tangent = 0.0;
#endif
}

float3 sample_aperiodic_tiling_far_albedo() {
  float3 face_mean = 0.0;

  [unroll]
  for (int face_id = 0; face_id < 10; ++face_id) {
    face_mean += aperiodic_face_color(face_id) * aperiodic_face_weights[face_id];
  }

  float inner_width = saturate(1.0 - 2.0 * min(_Aperiodic_Tiling_Edge_Thickness, 0.5));
  float inner_fraction = inner_width * inner_width;
  return lerp(_Aperiodic_Tiling_Edge_Color.rgb, face_mean, inner_fraction);
}

float aperiodic_minified_normal_kernel_roughness(float minified_weight) {
#if defined(_APERIODIC_TILING_NORMALS)
  float bevel_width = min(_Aperiodic_Tiling_Normal_Thickness, 0.5);
  if (bevel_width <= 1e-5 || minified_weight <= 1e-5) {
    return 0.0;
  }

  float flat_width = saturate(1.0 - 2.0 * bevel_width);
  float bevel_fraction = 1.0 - flat_width * flat_width;
  float slope_variance =
      _Aperiodic_Tiling_Normal_Strength * _Aperiodic_Tiling_Normal_Strength *
      bevel_fraction * APERIODIC_SMOOTHSTEP_SQ_MEAN;
  return minified_weight * min(2.0 * slope_variance, _Specular_AA_Threshold);
#else
  return 0.0;
#endif
}

void apply_aperiodic_smoothness(inout float smoothness, float minified_weight) {
  float kernel_roughness = aperiodic_minified_normal_kernel_roughness(minified_weight);
  if (kernel_roughness <= 1e-5) {
    return;
  }

  float perceptual_roughness = 1.0 - smoothness;
  float roughness = perceptual_roughness * perceptual_roughness;
  float square_roughness = saturate(roughness * roughness + kernel_roughness);
  smoothness = 1.0 - saturate(sqrt(sqrt(square_roughness)));
}

void sample_aperiodic_tiling_minified(float2 uv, float4 p03, float p44, float2 uv_ddx,
                                      float2 uv_ddy, float far_weight, out float3 albedo,
                                      out float3 tiling_normal_tangent) {
  float3 filtered_albedo;
  float3 filtered_normal_tangent;
  float3 far_albedo = sample_aperiodic_tiling_far_albedo();

  sample_aperiodic_tiling_filtered(
      uv, p03, p44, uv_ddx, uv_ddy, filtered_albedo, filtered_normal_tangent);
  albedo = lerp(filtered_albedo, far_albedo, far_weight);

#if defined(_APERIODIC_TILING_NORMALS)
  tiling_normal_tangent = normalize(lerp(
      filtered_normal_tangent, float3(0.0, 0.0, 1.0), far_weight));
#else
  tiling_normal_tangent = 0.0;
#endif
}

void sample_aperiodic_tiling(float2 uv, out float3 albedo, out float3 tiling_normal_tangent,
                             out float minified_weight) {
  float4 p03 = uv.x * basis_u5_03 + uv.y * basis_v5_03;
  float p44 = uv.x * basis_u5_44 + uv.y * basis_v5_44;
  float2 uv_ddx = ddx(uv);
  float2 uv_ddy = ddy(uv);
  AperiodicFootprintRange footprint = aperiodic_footprint_range(uv_ddx, uv_ddy);
  float filter_end = APERIODIC_FILTER_THRESHOLD + APERIODIC_FILTER_BLEND_WIDTH;
  float far_weight = smoothstep(
      APERIODIC_FILTER_THRESHOLD, APERIODIC_FAR_THRESHOLD, footprint.min_extent);

  [branch]
  if (footprint.max_extent >= filter_end) {
    sample_aperiodic_tiling_minified(
        uv, p03, p44, uv_ddx, uv_ddy, far_weight, albedo, tiling_normal_tangent);
    minified_weight = 1.0;
  } else {
    float3 point_albedo;
    float3 point_normal_tangent;
    float3 minified_albedo;
    float3 minified_normal_tangent;
    float filter_weight = smoothstep(0.0, filter_end, footprint.max_extent);

    sample_aperiodic_tiling_point(uv, p03, p44, point_albedo, point_normal_tangent);
    sample_aperiodic_tiling_minified(
        uv, p03, p44, uv_ddx, uv_ddy, far_weight, minified_albedo, minified_normal_tangent);
    albedo = lerp(point_albedo, minified_albedo, filter_weight);
    tiling_normal_tangent = normalize(lerp(
        point_normal_tangent, minified_normal_tangent, filter_weight));
    minified_weight = filter_weight;
  }
}
#endif  // defined(_APERIODIC_TILING)

void apply_aperiodic_tiling(float2 uv, inout float3 albedo, inout float smoothness,
                            inout float3 normal_tangent) {
#if defined(_APERIODIC_TILING)
  uv *= _Aperiodic_Tiling_Scale;
  float3 tiling_normal_tangent;
  float minified_weight;
  sample_aperiodic_tiling(uv, albedo, tiling_normal_tangent, minified_weight);
#if defined(_APERIODIC_TILING_NORMALS)
  normal_tangent = tiling_normal_tangent;
  apply_aperiodic_smoothness(smoothness, minified_weight);
#endif
#endif  // _APERIODIC_TILING
}

#endif  // __APERIODIC_TILING_INC