path: root/math.cginc

#ifndef __MATH_INC
#define __MATH_INC

#include "hashwithoutsine.cginc"

#define PI              3.14159265358979323846264f
#define TAU             (2.0f * PI)
#define HALF_PI         (PI * 0.5f)
#define RCP_PI          (1.0f / PI)
#define RCP_TAU         (1.0f / TAU)
#define PHI             1.618033989f
#define RCP_PHI         0.618033989f
#define SQRT_2          1.414213562f
#define SQRT_2_RCP      0.707106781f
#define RCP_SQRT_2      0.707106781f
#define RCP_SQRT_3      0.577350269f
#define TWO_OVER_THREE  0.6666666666666666f
#define TWO_OVER_SQRT_3 1.15470054f
#define SQRT_3          1.73205081f
#define SQRT_3_OVER_2   0.8660254037844386f
#define EULERS_CONSTANT 2.718281828f

#define F1_TO_F3(x) float3((x), (x), (x))

// Remaps [0, INT_MAX] to [0, 1]
#define UINT_TO_UNIT (1.0 / 4294967296.0)

float sin_noise_3d(float3 uvw) {
  return sin(uvw[0]) * sin(uvw[1]) * sin(uvw[2]);
}

float sin_noise_3d_fbm(float3 uvw, uint octaves, float k, float strength) {
  float result = 0;
  float factor = 1.0f;
  for (uint i = 0; i < octaves; i++) {
    result += sin_noise_3d(uvw) * factor * strength;
    uvw *= k;
    factor /= k;
  }
  return result;
}

// Wrap NoL. Assume it's already clamped.
// At k=0, you get standard lambertian shading.
// At k=0.5, you get half-lambertian shading.
// At k=1.0, you get flat shading.
// k must be on [0, 1].
// Energy preserving, within some small bound.
float wrapNoL(float NoL, float k) {
  float lambertian = NoL;
  float tmp = (NoL + 0.5f) / (1.0f + 0.5f);
  float half_lambertian = tmp * tmp;
  float flat = RCP_PI;

  if (k < 0.5) {
    return lerp(lambertian, half_lambertian, k * 2.0f);
  } else {
    return lerp(half_lambertian, flat, k * 2.0f - 1.0f);
  }
}

// Vector rotation with a quaternion
// https://blog.molecular-matters.com/2013/05/24/a-faster-quaternion-vector-multiplication/
float3 rotate_vector(float3 v, float4 q)
{
  float3 t = 2.0 * cross(q.xyz, v);
  return v + q.w * t + cross(q.xyz, t);
}

float4 get_quaternion(float3 axis_normal, float theta) {
  float sint2, cost2;
  sincos(theta*0.5, sint2, cost2);
  return float4(axis_normal * sint2, cost2);
}

// Cartesian to packed cube hexagonal coordinates.
// Stores (q + r, q, r), where the article's cube coordinates satisfy
// q + r + s = 0 and s = -(q + r).
// Based on this: https://backdrifting.net/post/064_hex_grids
float3 cart_to_hex(float2 cart) {
  float q = cart.x - cart.y * RCP_SQRT_3;
  float r = cart.y * (2.0f * RCP_SQRT_3);
  return float3(q + r, q, r);
}

float2 hex_to_cart(float3 hex_coord) {
  return float2(
      hex_coord[1] + hex_coord[2] * 0.5f,
      hex_coord[2] * SQRT_3_OVER_2);
}


float3 round_hex(float3 hex_coord) {
    float3 rounded = round(hex_coord);
    float3 diff = abs(rounded - hex_coord);
    float error = rounded.x - rounded.y - rounded.z;

    float3 max_mask = float3(
        diff.x > diff.y && diff.x > diff.z,
        diff.y > diff.x && diff.y > diff.z,
        diff.z > diff.x && diff.z > diff.y);

    rounded += error * float3(-1, 1, 1) * max_mask;
    return rounded;
}

// Reoriented normal mapping
// https://blog.selfshadow.com/publications/blending-in-detail/
// Inputs are in tangent space.
float3 blendNormalsHill12(float3 n0, float3 n1) {
  n0.z += 1.0;
  n1.xy = -n1.xy;
  return normalize(n0 * dot(n0, n1) - n1 * n0.z);
}

// https://en.wikipedia.org/wiki/Relative_luminance
float luminance(float3 rgb) {
  return dot(float3(0.2126, 0.7152, 0.0722), rgb);
}

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

// O'Neill-style PCG 32-bit output permutation (RXS-M-XS).
uint pcg32(uint input) {
  uint state = input * 747796405u + 2891336453u;
  uint word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
  return (word >> 22u) ^ word;
}

// 3 in 1 out hash. Mix three lanes in integer space, then PCG-finalize once.
uint hash31_u32(uint3 p) {
  uint seed = p.x * 0x2c1b3c6du;
  seed ^= p.y * 0x297a2d39u;
  seed ^= p.z * 0x1b56c4e9u;
  return pcg32(seed);
}

// Float wrapper for arbitrary inputs.
float hash31_ff(float3 p) {
  return hash31_u32(asuint(p)) * UINT_TO_UNIT;
}

float hash31_if(int3 p) {
  return hash31_u32(p) * UINT_TO_UNIT;
}

// Procedural value noise in [0,1]^3 — trilinear interpolation of hashed corners.
float3 value_noise3(float3 p) {
  int3 i = (int3)floor(p);
  float3 f = frac(p);
  float3 u = f * f * (3.0 - 2.0 * f);

  return lerp(
    lerp(lerp(hash31_if(i + int3(0, 0, 0)), hash31_if(i + int3(1, 0, 0)), u.x),
         lerp(hash31_if(i + int3(0, 1, 0)), hash31_if(i + int3(1, 1, 0)), u.x), u.y),
    lerp(lerp(hash31_if(i + int3(0, 0, 1)), hash31_if(i + int3(1, 0, 1)), u.x),
         lerp(hash31_if(i + int3(0, 1, 1)), hash31_if(i + int3(1, 1, 1)), u.x), u.y),
    u.z);
}

// Domain warping using procedural value noise.
float3 domain_warp_procedural(float3 uvw, float strength,
    uint octaves, float lacunarity, float gain) {
  float3 noise = 0;
  float g = 1;

  for (uint ii = 0; ii < octaves; ++ii) {
    // Need to remap noise onto [-1, 1] when domain warping.
    float3 vnoise = value_noise3(uvw + (noise * 2.0 - 1.0) * strength);
    noise += vnoise * g;
    uvw *= lacunarity;
    g *= gain;
  }

  noise *= (1 - gain) / (1 - pow(gain, octaves));

  return noise;
}

float3 value_noise_3d_tex(Texture3D tex, SamplerState s, float3 p) {
  float w, h, d;
  tex.GetDimensions(w, h, d);
  float3 res = float3(w, h, d);

  p *= res;
  float3 i = floor(p);
  float3 f = frac(p);
  float3 u = f * f * (3.0 - 2.0 * f);

  return tex.Sample(s, (i + 0.5 + u) / res).rgb;
}

// Domain warping using a 3D noise texture. Texture should have an EV of
// 0.5. Uses cubic interpolation between lattice points (same semantics as
// domain_warp_procedural / value_noise3).
float3 domain_warp_3d_tex(Texture3D noise_tex, SamplerState s, float3 uvw,
    float strength, uint octaves, float lacunarity, float gain) {
  float3 noise = 0;
  float g = 1;

  for (uint ii = 0; ii < octaves; ++ii) {
    noise += value_noise_3d_tex(noise_tex, s, uvw + noise * strength) * g;
    uvw *= lacunarity;
    g *= gain;
  }

  // Normalize: geometric series 1 + r + ... + r^{n-1} = (1 - r^n) / (1 - r)
  noise *= (1 - gain) / (1 - pow(gain, octaves));

  return noise;
}

// Return distance to the nearest voronoi cell edge. Also returns p1 and p2,
// vectors from the x to the nearest 2 lattice points, and p1_id, an identifier
// for p1 which does not within the current cell.
float voronoi_d_2d(float2 x, out float2 p1, out float2 p2, out float2 p1_id) {
  float2 x_floor = floor(x);
  float2 x_frac = frac(x);
  float d1 = 1e6;
  float d2 = 1e6;
  p1 = 0;
  p2 = 0;

  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      float2 cell_offset = float2(i, j);
      // Radial vector from origin lattice point to current.
      float2 r = cell_offset + hash22_fast(x_floor + cell_offset) - x_frac;
      float d = dot(r, r);
      if (d < d1) {
        d2 = d1;
        p2 = p1;
        d1 = d;
        p1 = r;
        p1_id = x_floor + cell_offset;
      } else if (d < d2) {
        d2 = d;
        p2 = r;
      }
    }
  }
  return (d2 - d1);
}

// stripe_dir is the direction we want to draw stripes along. It must be
// normalized.
float4 phacelle_noise_2d(float2 x, float2 stripe_dir) {
  float2 x_floor = floor(x);
  float2 x_frac = frac(x);

  float2 stripe_ortho = float2(-stripe_dir.y, stripe_dir.x) * TAU;

  float2 vec = 0;
  float normalization = 0;
  for (int j = -1; j <= 2; j++) {
    for (int i = -1; i <= 2; i++) {
      float2 cell_offset = float2(i, j);
      // Radial vector from origin lattice point to current.
      float2 r = cell_offset + (hash22_fast(x_floor + cell_offset) - 0.5) * 0.5 - x_frac;
      float d2 = dot(r, r);
      // 0th cell is centered at 0.5.
      float weight = max(0, exp(-2.0 * d2) - 0.01111);
      normalization += weight;
      float wave_input = dot(r, stripe_ortho);
      // TODO what the fuck
      // why is r not normalized
      float vecs, vecc;
      sincos(wave_input, vecs, vecc);
      vec += float2(vecc, vecs) * weight;
    }
  }
  vec /= normalization;
  float magnitude = dot(vec, vec);
  return float4(vec / magnitude, stripe_ortho);
}

float voronoi_d_2d(float2 x) {
  float2 p1, p2, p1_id;
  return voronoi_d_2d(x, p1, p2, p1_id);
}

// Return distance to the nearest voronoi cell edge, clamped to [0, 0.5].
// 0.5 is on the edge, 0 is far from it.
float voronoi_d_3d(float3 x) {
  float3 x_floor = floor(x);
  float3 x_frac = frac(x);
  float d1 = 1e6;
  float d2 = 1e6;
  float3 p1 = 0;
  float3 p2 = 0;

  for (int k = -1; k <= 1; k++) {
    for (int j = -1; j <= 1; j++) {
      for (int i = -1; i <= 1; i++) {
        float3 cell_offset = float3(i, j, k);
        float3 r = cell_offset + hash33_fast(x_floor + cell_offset) - x_frac;
        float d = dot(r, r);
        if (d < d1) {
          d2 = d1;
          p2 = p1;
          d1 = d;
          p1 = r;
        } else if (d < d2) {
          d2 = d;
          p2 = r;
        }
      }
    }
  }
  float d = (d2 - d1) / (2.0f * max(1e-4, length(p2 - p1)));
  return max(0.0f, 0.5f - d);
}

float median(float3 x) {
  // Get the min and max.
  float x_min= min(min(x.r, x.g), x.b);
  float x_max = max(max(x.r, x.g), x.b);

  // Compute (x.r + x.g + x.b) - (x_min + x_max). This gives us the median.
  return (x.r + x.g + x.b) - (x_min + x_max);
}

float3 linear_to_srgb(float3 linear_color) {
  float3 lo = 12.92f * linear_color;
  float3 hi = 1.055f * pow(linear_color, 1.0f / 2.4f) - 0.055f;
  return lerp(lo, hi, step(0.0031308f, linear_color));
}

float3 srgb_to_linear(float3 srgb_color) {
  float3 lo = srgb_color / 12.92f;
  float3 hi = pow((srgb_color + 0.055f) / 1.055f, 2.4f);
  return lerp(lo, hi, step(0.04045f, srgb_color));
}

#endif  // __MATH_INC