Finish implementing burley per-channel histogram preserving blending operator

1 files changed, 112 insertions, 15 deletions

diff --git a/pbr.cginc b/pbr.cginc
index d59001f..ecfd083 100755
--- a/pbr.cginc
+++ b/pbr.cginc
@@ -181,34 +181,131 @@ void apply_letter_grid(v2f i, inout Pbr pbr) {
 #endif
 }
 
-void apply_burley_tiling(v2f i, inout Pbr pbr) {
 #if defined(_BURLEY_TILING)
-  float2 uv = i.uv01.xy - 0.5;
-  uv *= TWO_OVER_SQRT_3;
-  uv /= _Burley_Tiling_Output_Scale;
-  float3 hex_coord = cart_to_hex(uv);
-  float3 cell = round_hex(hex_coord);
-  float3 cube_id = float3(cell.y, cell.z, -cell.y - cell.z);
-  float3 local_hex = hex_coord - cell;
-  // Get UVs on [-1/2,1/2] for the current cell.
-  // This is the tightest square fit that preserves the hex aspect ratio.
-  float2 local_uv = hex_to_cart(local_hex) * SQRT_3_OVER_2;
+float2 burley_tri_to_cart(float2 tri_coord) {
+  return float2(
+      tri_coord.x + tri_coord.y * 0.5f,
+      tri_coord.y * SQRT_3_OVER_2);
+}
+
+float3 burley_apply_blend_gamma(float3 weights, float gamma) {
+  weights = pow(weights, gamma);
+  return weights / (weights.x + weights.y + weights.z);
+}
+
+// Equation 4 (first half).
+float3 burley_soft_clipping_lower_half(float3 x_hat, float w_hat) {
+  float linear_start = 0.25f * (2.0f - w_hat);
+  float3 linear_value = (x_hat - 0.5f) / w_hat + 0.5f;
+  float3 linear_mask = step(float3(linear_start, linear_start, linear_start), x_hat);
+
+  if (w_hat >= TWO_OVER_THREE) {
+    float3 t = x_hat / (2.0f - w_hat);
+    float3 quadratic = 8.0f * (1.0f / w_hat - 1.0f) * t * t + (3.0f - 2.0f / w_hat) * t;
+    return lerp(quadratic, linear_value, linear_mask);
+  }
+
+  float quadratic_start = 0.25f * (2.0f - 3.0f * w_hat);
+  float3 d = (x_hat - quadratic_start) / w_hat;
+  float3 quadratic = d * d;
+  float3 quadratic_mask = step(float3(quadratic_start, quadratic_start, quadratic_start), x_hat);
+  float3 result = quadratic * quadratic_mask;
+  return lerp(result, linear_value, linear_mask);
+}
+
+// Equation 4.
+float3 burley_soft_clipping_contrast(float3 x_hat, float w_hat) {
+  float3 upper_mask = step(0.5f, x_hat);
+  float3 lower_x = min(x_hat, 1.0f - x_hat);
+  float3 lower_y = burley_soft_clipping_lower_half(lower_x, w_hat);
+  return lerp(lower_y, 1.0f - lower_y, upper_mask);
+}
+
+float3 burley_apply_soft_clipping(float3 gaussian_color, float3 weights) {
+  float w_hat = sqrt(dot(weights, weights));
+  return burley_soft_clipping_contrast(gaussian_color, w_hat);
+}
+
+float3 burley_degaussianize(float3 gaussian_color) {
+  float2 uv_r = float2(gaussian_color.r, 0.5f);
+  float2 uv_g = float2(gaussian_color.g, 0.5f);
+  float2 uv_b = float2(gaussian_color.b, 0.5f);
+  float lut_r = _Burley_Tiling_Maintex_LUT.Sample(linear_clamp_s, uv_r).r;
+  float lut_g = _Burley_Tiling_Maintex_LUT.Sample(linear_clamp_s, uv_g).g;
+  float lut_b = _Burley_Tiling_Maintex_LUT.Sample(linear_clamp_s, uv_b).b;
+  return srgb_to_linear(float3(lut_r, lut_g, lut_b));
+}
+
+float4 burley_sample_patch(float2 uv, float2 uv_dx, float2 uv_dy, float2 tri_vertex,
+                           float input_scale) {
+  float3 cube_id = float3(tri_vertex.x, tri_vertex.y, -tri_vertex.x - tri_vertex.y);
   float3 tile_rand3 = hash33_fast(cube_id);
+  float2 vertex_uv = burley_tri_to_cart(tri_vertex);
+  // Map the unit-radius hex support to the unit square so arbitrary rotation
+  // stays within bounds.
+  float2 local_uv = (uv - vertex_uv) * 0.5f;
   // Apply input scaling.
-  float input_scale = saturate(_Burley_Tiling_Input_Scale);
   local_uv *= input_scale;
+  float2 sample_dx = uv_dx * (0.5f * input_scale);
+  float2 sample_dy = uv_dy * (0.5f * input_scale);
   // Rotate.
   float theta = hash31_ff(tile_rand3) * TAU;
   float2x2 rot = float2x2(cos(theta), -sin(theta), sin(theta), cos(theta));
   local_uv = mul(rot, local_uv);
+  sample_dx = mul(rot, sample_dx);
+  sample_dy = mul(rot, sample_dy);
   // Apply randomized offset, staying within bounds.
   // The scaled-and-rotated footprint is bounded by [-Input_Scale / 2, Input_Scale / 2],
   // so we can offset by [(1 - Input_Scale) / 2].
-  float2 random_offset = (tile_rand3.yz * 2.0 - 1.0) * (0.5 * (1.0 - input_scale));
+  float2 random_offset = (tile_rand3.yz * 2.0f - 1.0f) * (0.5f * (1.0f - input_scale));
   local_uv += random_offset;
   // Finally, remap onto [0, 1].
-  local_uv += 0.5;
-  pbr.albedo.xy = local_uv.xy;
+  local_uv += 0.5f;
+
+  float2 sample_uv = local_uv * _Burley_Tiling_Maintex_ST.xy + _Burley_Tiling_Maintex_ST.zw;
+  sample_dx *= _Burley_Tiling_Maintex_ST.xy;
+  sample_dy *= _Burley_Tiling_Maintex_ST.xy;
+  return _Burley_Tiling_Maintex.SampleGrad(
+      aniso4_trilinear_repeat_s, sample_uv, sample_dx, sample_dy);
+}
+#endif  // _BURLEY_TILING
+
+void apply_burley_tiling(v2f i, inout Pbr pbr) {
+#if defined(_BURLEY_TILING)
+  // Center at 0.
+  float2 uv = i.uv01.xy - 0.5;
+  // Scale so that any rotation remains within [0, 1] bounds.
+  uv *= TWO_OVER_SQRT_3;
+  uv /= _Burley_Tiling_Output_Scale;
+  float3 hex_coord = cart_to_hex(uv);
+  float2 tri_coord = hex_coord.yz;
+  float2 tri_cell = floor(tri_coord);
+  float2 tri_frac = tri_coord - tri_cell;
+  float2 vertex_0;
+  float2 vertex_1;
+  float2 vertex_2;
+  float3 baryc;
+  if (tri_frac.x + tri_frac.y < 1.0f) {
+    vertex_0 = tri_cell;
+    vertex_1 = tri_cell + float2(1.0f, 0.0f);
+    vertex_2 = tri_cell + float2(0.0f, 1.0f);
+    baryc = float3(1.0f - (tri_frac.x + tri_frac.y), tri_frac.x, tri_frac.y);
+  } else {
+    vertex_0 = tri_cell + 1.0f;
+    vertex_1 = tri_cell + float2(0.0f, 1.0f);
+    vertex_2 = tri_cell + float2(1.0f, 0.0f);
+    baryc = float3(tri_frac.x + tri_frac.y - 1.0f, 1.0f - tri_frac.x, 1.0f - tri_frac.y);
+  }
+
+  float input_scale = _Burley_Tiling_Input_Scale;
+  float3 weights = burley_apply_blend_gamma(baryc, _Burley_Tiling_Blend_Gamma);
+  float2 uv_dx = ddx(uv);
+  float2 uv_dy = ddy(uv);
+  float4 patch_0 = burley_sample_patch(uv, uv_dx, uv_dy, vertex_0, input_scale);
+  float4 patch_1 = burley_sample_patch(uv, uv_dx, uv_dy, vertex_1, input_scale);
+  float4 patch_2 = burley_sample_patch(uv, uv_dx, uv_dy, vertex_2, input_scale);
+  float4 gaussian_blend = patch_0 * weights.x + patch_1 * weights.y + patch_2 * weights.z;
+  pbr.albedo.xyz = burley_degaussianize(burley_apply_soft_clipping(gaussian_blend.rgb, weights));
 #endif  // _BURLEY_TILING
 }