e89fd68778
The cyclic non-strict step(x0.yzx, x0.xyz) made the four corners a non-simplex where x0's components are equal: the origin and every (t, t, t) read a different function (snoise(0) was -0.4358730, jumps up to 0.77). One strict compare, z against x, gives simplexnoise1234's order. Off the tie set old and new return the same bits. NoiseGradTests.SnoiseIsContinuousOnTheSimplexDiagonal fails against the old rule (0.054 within 1 ulp of (-53.6, -53.6, -53.6)). The origin golden is 0 and a (1.25, 1.25, 1.25) diagonal row joins it. SNOISE_LIPSCHITZ and the GLSL follow in the next commits, so the ratio pin and CommittedGlslEqualsLowering are red here. Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
1023 lines
56 KiB
C#
1023 lines
56 KiB
C#
using System;
|
|
|
|
// Cut 2a-i fix batch 1: cellular_unit's [0,1) mapping bug (F7) is only reachable from cellular's
|
|
// private call sites; this lets CultMath.Tests probe it directly instead of hunting for a hash
|
|
// input that happens to land on the bit pattern under test.
|
|
[assembly: System.Runtime.CompilerServices.InternalsVisibleTo("CultMath.Tests")]
|
|
|
|
namespace CultMath;
|
|
|
|
public static partial class math
|
|
{
|
|
public const float PI = MathF.PI;
|
|
public const float TAU = MathF.PI * 2.0f;
|
|
public const float HALF_PI = MathF.PI * 0.5f;
|
|
|
|
/// <summary>
|
|
/// Extrinsic Euler order for <see cref="CultMath.float3x3.Euler(CultMath.float3, RotationOrder)"/>:
|
|
/// the first letter is the axis applied first. Nested here so <c>using static CultMath.math;</c>
|
|
/// brings it into scope with the rest of the shader vocabulary.
|
|
/// </summary>
|
|
public enum RotationOrder : byte
|
|
{
|
|
XYZ,
|
|
XZY,
|
|
YXZ,
|
|
YZX,
|
|
ZXY,
|
|
ZYX,
|
|
Default = ZXY,
|
|
}
|
|
|
|
public static float2 float2(float x, float y) => new(x, y);
|
|
public static float2 float2(float value) => new(value, value);
|
|
public static bool2 bool2(bool x, bool y) => new(x, y);
|
|
public static bool2 bool2(bool value) => new(value, value);
|
|
public static int2 int2(int x, int y) => new(x, y);
|
|
public static int2 int2(int value) => new(value, value);
|
|
// Mixed constructor functions (float4(float2, float2), ...) are generated in Swizzles.g.cs.
|
|
public static float3 float3(float x, float y, float z) => new(x, y, z);
|
|
public static float3 float3(float value) => new(value, value, value);
|
|
public static float4 float4(float x, float y, float z, float w) => new(x, y, z, w);
|
|
public static float4 float4(float value) => new(value, value, value, value);
|
|
public static bool3 bool3(bool x, bool y, bool z) => new(x, y, z);
|
|
public static bool3 bool3(bool value) => new(value, value, value);
|
|
public static bool4 bool4(bool x, bool y, bool z, bool w) => new(x, y, z, w);
|
|
public static bool4 bool4(bool value) => new(value, value, value, value);
|
|
public static int3 int3(int x, int y, int z) => new(x, y, z);
|
|
public static int3 int3(int value) => new(value, value, value);
|
|
public static int4 int4(int x, int y, int z, int w) => new(x, y, z, w);
|
|
public static int4 int4(int value) => new(value, value, value, value);
|
|
|
|
// Matrix constructor functions take rows, as HLSL does. Under `using static CultMath.math;`
|
|
// these method groups hide the type names in member access, so static members need a
|
|
// qualified type: `CultMath.float3x3.Euler(...)`, `CultMath.float2x2.Rotate(...)`.
|
|
public static float2x2 float2x2(float m00, float m01, float m10, float m11) => new(m00, m01, m10, m11);
|
|
public static float2x2 float2x2(float2 row0, float2 row1) => new(row0, row1);
|
|
public static float3x3 float3x3(float m00, float m01, float m02, float m10, float m11, float m12, float m20, float m21, float m22) =>
|
|
new(m00, m01, m02, m10, m11, m12, m20, m21, m22);
|
|
public static float3x3 float3x3(float3 row0, float3 row1, float3 row2) => new(row0, row1, row2);
|
|
public static float radians(float degrees) => degrees * (PI / 180.0f);
|
|
public static float degrees(float radians) => radians * (180.0f / PI);
|
|
public static float sin(float value) => MathF.Sin(value);
|
|
public static float cos(float value) => MathF.Cos(value);
|
|
public static float tan(float value) => MathF.Tan(value);
|
|
public static float asin(float value) => MathF.Asin(value);
|
|
public static float acos(float value) => MathF.Acos(value);
|
|
public static float atan(float value) => MathF.Atan(value);
|
|
public static float atan2(float y, float x) => MathF.Atan2(y, x);
|
|
public static float sqrt(float value) => MathF.Sqrt(value);
|
|
public static float pow(float value, float power) => MathF.Pow(value, power);
|
|
public static bool isinf(float value) => float.IsInfinity(value);
|
|
public static float exp(float value) => MathF.Exp(value);
|
|
public static double exp(double value) => Math.Exp(value);
|
|
public static float2 exp(float2 value) => new(exp(value.x), exp(value.y));
|
|
public static float3 exp(float3 value) => new(exp(value.x), exp(value.y), exp(value.z));
|
|
public static float4 exp(float4 value) => new(exp(value.x), exp(value.y), exp(value.z), exp(value.w));
|
|
public static double2 exp(double2 value) => new(exp(value.x), exp(value.y));
|
|
public static double3 exp(double3 value) => new(exp(value.x), exp(value.y), exp(value.z));
|
|
|
|
public static float abs(float value) => MathF.Abs(value);
|
|
public static float2 abs(float2 value) => new(abs(value.x), abs(value.y));
|
|
public static float3 abs(float3 value) => new(abs(value.x), abs(value.y), abs(value.z));
|
|
public static float4 abs(float4 value) => new(abs(value.x), abs(value.y), abs(value.z), abs(value.w));
|
|
public static int abs(int value) => Math.Abs(value);
|
|
public static int2 abs(int2 value) => new(abs(value.x), abs(value.y));
|
|
public static int3 abs(int3 value) => new(abs(value.x), abs(value.y), abs(value.z));
|
|
public static int4 abs(int4 value) => new(abs(value.x), abs(value.y), abs(value.z), abs(value.w));
|
|
|
|
// HLSL sign returns int. NaN compares false both ways, so it yields 0 instead of throwing.
|
|
public static int sign(float value) => value > 0.0f ? 1 : value < 0.0f ? -1 : 0;
|
|
public static int2 sign(float2 value) => new(sign(value.x), sign(value.y));
|
|
public static int3 sign(float3 value) => new(sign(value.x), sign(value.y), sign(value.z));
|
|
public static int4 sign(float4 value) => new(sign(value.x), sign(value.y), sign(value.z), sign(value.w));
|
|
public static int sign(int value) => value > 0 ? 1 : value < 0 ? -1 : 0;
|
|
public static int2 sign(int2 value) => new(sign(value.x), sign(value.y));
|
|
public static int3 sign(int3 value) => new(sign(value.x), sign(value.y), sign(value.z));
|
|
public static int4 sign(int4 value) => new(sign(value.x), sign(value.y), sign(value.z), sign(value.w));
|
|
|
|
public static float floor(float value) => MathF.Floor(value);
|
|
public static float2 floor(float2 value) => new(floor(value.x), floor(value.y));
|
|
public static float3 floor(float3 value) => new(floor(value.x), floor(value.y), floor(value.z));
|
|
public static float4 floor(float4 value) => new(floor(value.x), floor(value.y), floor(value.z), floor(value.w));
|
|
|
|
public static float ceil(float value) => MathF.Ceiling(value);
|
|
public static float2 ceil(float2 value) => new(ceil(value.x), ceil(value.y));
|
|
public static float3 ceil(float3 value) => new(ceil(value.x), ceil(value.y), ceil(value.z));
|
|
public static float4 ceil(float4 value) => new(ceil(value.x), ceil(value.y), ceil(value.z), ceil(value.w));
|
|
|
|
public static float frac(float value) => value - floor(value);
|
|
public static float2 frac(float2 value) => value - floor(value);
|
|
public static float3 frac(float3 value) => value - floor(value);
|
|
public static float4 frac(float4 value) => value - floor(value);
|
|
// dxc has no double Frc: it demotes to float, as it does for exp and lerp. CultMath's double
|
|
// overloads keep double precision instead, so frac(double) is x - floor(x) in double.
|
|
public static double frac(double value) => value - Math.Floor(value);
|
|
|
|
// dxc lowers float min/max/clamp to DXIL FMin/FMax (saturate to Saturate, defined as FMin(1, FMax(0, x))).
|
|
// DXIL.rst: FMin(a, b) is a < b ? a : b, FMax(a, b) is a >= b ? a : b, and a NaN operand returns the other.
|
|
public static float min(float left, float right) => left < right || float.IsNaN(right) ? left : right;
|
|
public static int min(int left, int right) => Math.Min(left, right);
|
|
public static float2 min(float2 left, float2 right) => new(min(left.x, right.x), min(left.y, right.y));
|
|
public static float3 min(float3 left, float3 right) => new(min(left.x, right.x), min(left.y, right.y), min(left.z, right.z));
|
|
public static float4 min(float4 left, float4 right) => new(min(left.x, right.x), min(left.y, right.y), min(left.z, right.z), min(left.w, right.w));
|
|
public static int2 min(int2 left, int2 right) => new(min(left.x, right.x), min(left.y, right.y));
|
|
public static int3 min(int3 left, int3 right) => new(min(left.x, right.x), min(left.y, right.y), min(left.z, right.z));
|
|
public static int4 min(int4 left, int4 right) => new(min(left.x, right.x), min(left.y, right.y), min(left.z, right.z), min(left.w, right.w));
|
|
|
|
public static float max(float left, float right) => left >= right || float.IsNaN(right) ? left : right;
|
|
public static int max(int left, int right) => Math.Max(left, right);
|
|
public static float2 max(float2 left, float2 right) => new(max(left.x, right.x), max(left.y, right.y));
|
|
public static float3 max(float3 left, float3 right) => new(max(left.x, right.x), max(left.y, right.y), max(left.z, right.z));
|
|
public static float4 max(float4 left, float4 right) => new(max(left.x, right.x), max(left.y, right.y), max(left.z, right.z), max(left.w, right.w));
|
|
public static int2 max(int2 left, int2 right) => new(max(left.x, right.x), max(left.y, right.y));
|
|
public static int3 max(int3 left, int3 right) => new(max(left.x, right.x), max(left.y, right.y), max(left.z, right.z));
|
|
public static int4 max(int4 left, int4 right) => new(max(left.x, right.x), max(left.y, right.y), max(left.z, right.z), max(left.w, right.w));
|
|
|
|
public static float clamp(float value, float minimum, float maximum) => min(max(value, minimum), maximum);
|
|
public static float2 clamp(float2 value, float2 minimum, float2 maximum) => min(max(value, minimum), maximum);
|
|
public static float3 clamp(float3 value, float3 minimum, float3 maximum) => min(max(value, minimum), maximum);
|
|
public static float4 clamp(float4 value, float4 minimum, float4 maximum) => min(max(value, minimum), maximum);
|
|
// dxc lowers int clamp to IMin(IMax(x, a), b): inverted bounds return b rather than throwing like Math.Clamp.
|
|
public static int clamp(int value, int minimum, int maximum) => min(max(value, minimum), maximum);
|
|
public static int2 clamp(int2 value, int2 minimum, int2 maximum) => min(max(value, minimum), maximum);
|
|
public static int3 clamp(int3 value, int3 minimum, int3 maximum) => min(max(value, minimum), maximum);
|
|
public static int4 clamp(int4 value, int4 minimum, int4 maximum) => min(max(value, minimum), maximum);
|
|
|
|
// DXIL Saturate is FMin(1, FMax(0, x)); that operand order, unlike clamp(x, 0, 1), maps -0 to +0.
|
|
public static float saturate(float value) => min(1.0f, max(0.0f, value));
|
|
public static float2 saturate(float2 value) => min(1.0f, max(0.0f, value));
|
|
public static float3 saturate(float3 value) => min(1.0f, max(0.0f, value));
|
|
public static float4 saturate(float4 value) => min(1.0f, max(0.0f, value));
|
|
|
|
public static float lerp(float start, float end, float amount) => start + (end - start) * amount;
|
|
public static double lerp(double start, double end, double amount) => start + (end - start) * amount;
|
|
public static float2 lerp(float2 start, float2 end, float2 amount) => start + (end - start) * amount;
|
|
public static float3 lerp(float3 start, float3 end, float3 amount) => start + (end - start) * amount;
|
|
public static float4 lerp(float4 start, float4 end, float4 amount) => start + (end - start) * amount;
|
|
public static double2 lerp(double2 start, double2 end, double2 amount) => start + (end - start) * amount;
|
|
public static double3 lerp(double3 start, double3 end, double3 amount) => start + (end - start) * amount;
|
|
|
|
public static float unlerp(float start, float end, float value) => (value - start) / (end - start);
|
|
public static double unlerp(double start, double end, double value) => (value - start) / (end - start);
|
|
public static float2 unlerp(float2 start, float2 end, float2 value) => (value - start) / (end - start);
|
|
public static float3 unlerp(float3 start, float3 end, float3 value) => (value - start) / (end - start);
|
|
public static float4 unlerp(float4 start, float4 end, float4 value) => (value - start) / (end - start);
|
|
public static double2 unlerp(double2 start, double2 end, double2 value) => (value - start) / (end - start);
|
|
public static double3 unlerp(double3 start, double3 end, double3 value) => (value - start) / (end - start);
|
|
|
|
// dxc lowers step(y, x) to x < y ? 0 : 1 (fcmp olt, then select), so a NaN on either side yields 1.
|
|
public static float step(float edge, float value) => value < edge ? 0.0f : 1.0f;
|
|
public static float2 step(float2 edge, float2 value) => new(step(edge.x, value.x), step(edge.y, value.y));
|
|
public static float3 step(float3 edge, float3 value) => new(step(edge.x, value.x), step(edge.y, value.y), step(edge.z, value.z));
|
|
public static float4 step(float4 edge, float4 value) => new(step(edge.x, value.x), step(edge.y, value.y), step(edge.z, value.z), step(edge.w, value.w));
|
|
|
|
public static float smoothstep(float minimum, float maximum, float value)
|
|
{
|
|
var t = saturate((value - minimum) / (maximum - minimum));
|
|
return t * t * (3.0f - 2.0f * t);
|
|
}
|
|
|
|
public static float2 smoothstep(float2 minimum, float2 maximum, float2 value)
|
|
{
|
|
var t = saturate((value - minimum) / (maximum - minimum));
|
|
return t * t * (3.0f - 2.0f * t);
|
|
}
|
|
|
|
public static float3 smoothstep(float3 minimum, float3 maximum, float3 value)
|
|
{
|
|
var t = saturate((value - minimum) / (maximum - minimum));
|
|
return t * t * (3.0f - 2.0f * t);
|
|
}
|
|
|
|
public static float4 smoothstep(float4 minimum, float4 maximum, float4 value)
|
|
{
|
|
var t = saturate((value - minimum) / (maximum - minimum));
|
|
return t * t * (3.0f - 2.0f * t);
|
|
}
|
|
|
|
public static float dot(float2 left, float2 right) => left.x * right.x + left.y * right.y;
|
|
public static float dot(float3 left, float3 right) => left.x * right.x + left.y * right.y + left.z * right.z;
|
|
public static float dot(float4 left, float4 right) => left.x * right.x + left.y * right.y + left.z * right.z + left.w * right.w;
|
|
|
|
public static float3 cross(float3 left, float3 right) =>
|
|
new(
|
|
left.y * right.z - left.z * right.y,
|
|
left.z * right.x - left.x * right.z,
|
|
left.x * right.y - left.y * right.x);
|
|
|
|
public static float lengthsq(float2 value) => dot(value, value);
|
|
public static float lengthsq(float3 value) => dot(value, value);
|
|
public static float lengthsq(float4 value) => dot(value, value);
|
|
|
|
public static float length(float2 value) => MathF.Sqrt(lengthsq(value));
|
|
public static float length(float3 value) => MathF.Sqrt(lengthsq(value));
|
|
public static float length(float4 value) => MathF.Sqrt(lengthsq(value));
|
|
|
|
public static float distance(float2 left, float2 right) => length(left - right);
|
|
public static float distance(float3 left, float3 right) => length(left - right);
|
|
public static float distance(float4 left, float4 right) => length(left - right);
|
|
|
|
// dxc lowers normalize(x) to x * Rsqrt(Dot(x, x)), and DXIL Rsqrt is 1 / sqrt(src): a zero vector is 0 * inf = NaN.
|
|
public static float2 normalize(float2 value) => value * (1.0f / MathF.Sqrt(dot(value, value)));
|
|
public static float3 normalize(float3 value) => value * (1.0f / MathF.Sqrt(dot(value, value)));
|
|
public static float4 normalize(float4 value) => value * (1.0f / MathF.Sqrt(dot(value, value)));
|
|
public static quaternion normalize(quaternion value)
|
|
{
|
|
var length = MathF.Sqrt(value.x * value.x + value.y * value.y + value.z * value.z + value.w * value.w);
|
|
if (length <= 1.0e-20f)
|
|
{
|
|
return quaternion.identity;
|
|
}
|
|
|
|
var invLength = 1.0f / length;
|
|
return new quaternion(value.x * invLength, value.y * invLength, value.z * invLength, value.w * invLength);
|
|
}
|
|
|
|
public static float2 reflect(float2 incident, float2 normal) => incident - 2.0f * dot(normal, incident) * normal;
|
|
public static float3 reflect(float3 incident, float3 normal) => incident - 2.0f * dot(normal, incident) * normal;
|
|
|
|
public static float2 rotate(float2 value, float radians)
|
|
{
|
|
var s = sin(radians);
|
|
var c = cos(radians);
|
|
return new float2(value.x * c - value.y * s, value.x * s + value.y * c);
|
|
}
|
|
|
|
public static float2 rotate_degrees(float2 value, float degrees) => rotate(value, radians(degrees));
|
|
|
|
public static float csum(float2 value) => value.x + value.y;
|
|
public static float csum(float3 value) => value.x + value.y + value.z;
|
|
public static float csum(float4 value) => value.x + value.y + value.z + value.w;
|
|
|
|
public static float decay(float source, float lambda, float dt) => source * exp(-lambda * dt);
|
|
public static float2 decay(float2 source, float lambda, float dt) => source * exp(-lambda * dt);
|
|
public static float3 decay(float3 source, float lambda, float dt) => source * exp(-lambda * dt);
|
|
public static float4 decay(float4 source, float lambda, float dt) => source * exp(-lambda * dt);
|
|
public static float2 decay(float2 source, float2 lambda, float dt) => source * exp(-lambda * dt);
|
|
public static float3 decay(float3 source, float3 lambda, float dt) => source * exp(-lambda * dt);
|
|
public static float4 decay(float4 source, float4 lambda, float dt) => source * exp(-lambda * dt);
|
|
|
|
public static float damp(float start, float end, float lambda, float dt) => lerp(start, end, 1.0f - exp(-lambda * dt));
|
|
public static float2 damp(float2 start, float2 end, float lambda, float dt) => lerp(start, end, 1.0f - exp(-lambda * dt));
|
|
public static float3 damp(float3 start, float3 end, float lambda, float dt) => lerp(start, end, 1.0f - exp(-lambda * dt));
|
|
public static float4 damp(float4 start, float4 end, float lambda, float dt) => lerp(start, end, 1.0f - exp(-lambda * dt));
|
|
|
|
public static float first_order_intercept_time(float shotSpeed, float3 targetRelativePosition, float3 targetRelativeVelocity)
|
|
{
|
|
var velocitySquared = lengthsq(targetRelativeVelocity);
|
|
if (velocitySquared < 0.001f)
|
|
{
|
|
return 0.0f;
|
|
}
|
|
|
|
var a = velocitySquared - shotSpeed * shotSpeed;
|
|
if (abs(a) < 0.001f)
|
|
{
|
|
var time = -lengthsq(targetRelativePosition) / (2.0f * dot(targetRelativeVelocity, targetRelativePosition));
|
|
return max(time, 0.0f);
|
|
}
|
|
|
|
var b = 2.0f * dot(targetRelativeVelocity, targetRelativePosition);
|
|
var c = lengthsq(targetRelativePosition);
|
|
var determinant = b * b - 4.0f * a * c;
|
|
|
|
if (determinant > 0.0f)
|
|
{
|
|
var root = sqrt(determinant);
|
|
var t1 = (-b + root) / (2.0f * a);
|
|
var t2 = (-b - root) / (2.0f * a);
|
|
if (t1 > 0.0f)
|
|
{
|
|
return t2 > 0.0f ? min(t1, t2) : t1;
|
|
}
|
|
|
|
return max(t2, 0.0f);
|
|
}
|
|
|
|
if (determinant < 0.0f)
|
|
{
|
|
return 0.0f;
|
|
}
|
|
|
|
return max(-b / (2.0f * a), 0.0f);
|
|
}
|
|
|
|
public static float3 first_order_intercept(
|
|
float3 shooterPosition,
|
|
float3 shooterVelocity,
|
|
float shotSpeed,
|
|
float3 targetPosition,
|
|
float3 targetVelocity)
|
|
{
|
|
var targetRelativePosition = targetPosition - shooterPosition;
|
|
var targetRelativeVelocity = targetVelocity - shooterVelocity;
|
|
var time = first_order_intercept_time(shotSpeed, targetRelativePosition, targetRelativeVelocity);
|
|
return targetPosition + time * targetRelativeVelocity;
|
|
}
|
|
|
|
public static float distance_to_segment(float2 point, float2 start, float2 end, out float2 closest)
|
|
{
|
|
var segment = end - start;
|
|
var segmentLengthSq = lengthsq(segment);
|
|
if (segmentLengthSq <= 0.0f)
|
|
{
|
|
closest = start;
|
|
return distance(point, start);
|
|
}
|
|
|
|
var t = saturate(dot(point - start, segment) / segmentLengthSq);
|
|
closest = start + segment * t;
|
|
return distance(point, closest);
|
|
}
|
|
|
|
public static float catmullrom(float p0, float p1, float p2, float p3, float t)
|
|
{
|
|
var t2 = t * t;
|
|
var t3 = t2 * t;
|
|
return 0.5f * ((2.0f * p1) + (p2 - p0) * t + (2.0f * p0 - 5.0f * p1 + 4.0f * p2 - p3) * t2 + (-p0 + 3.0f * p1 - 3.0f * p2 + p3) * t3);
|
|
}
|
|
|
|
public static float2 catmullrom(float2 p0, float2 p1, float2 p2, float2 p3, float t) =>
|
|
0.5f * ((2.0f * p1) + (p2 - p0) * t + (2.0f * p0 - 5.0f * p1 + 4.0f * p2 - p3) * (t * t) + (-p0 + 3.0f * p1 - 3.0f * p2 + p3) * (t * t * t));
|
|
|
|
public static float3 catmullrom(float3 p0, float3 p1, float3 p2, float3 p3, float t) =>
|
|
0.5f * ((2.0f * p1) + (p2 - p0) * t + (2.0f * p0 - 5.0f * p1 + 4.0f * p2 - p3) * (t * t) + (-p0 + 3.0f * p1 - 3.0f * p2 + p3) * (t * t * t));
|
|
|
|
public static float4 catmullrom(float4 p0, float4 p1, float4 p2, float4 p3, float t) =>
|
|
0.5f * ((2.0f * p1) + (p2 - p0) * t + (2.0f * p0 - 5.0f * p1 + 4.0f * p2 - p3) * (t * t) + (-p0 + 3.0f * p1 - 3.0f * p2 + p3) * (t * t * t));
|
|
|
|
public static float quadratic_bezier(float p0, float p1, float p2, float t)
|
|
{
|
|
t = saturate(t);
|
|
var inv = 1.0f - t;
|
|
return inv * inv * p0 + 2.0f * inv * t * p1 + t * t * p2;
|
|
}
|
|
|
|
public static float3 quadratic_bezier(float3 p0, float3 p1, float3 p2, float t)
|
|
{
|
|
t = saturate(t);
|
|
var inv = 1.0f - t;
|
|
return inv * inv * p0 + 2.0f * inv * t * p1 + t * t * p2;
|
|
}
|
|
|
|
public static float cubic_bezier(float p0, float p1, float p2, float p3, float t)
|
|
{
|
|
t = saturate(t);
|
|
var inv = 1.0f - t;
|
|
return inv * inv * inv * p0 + 3.0f * inv * inv * t * p1 + 3.0f * inv * t * t * p2 + t * t * t * p3;
|
|
}
|
|
|
|
public static float3 cubic_bezier(float3 p0, float3 p1, float3 p2, float3 p3, float t)
|
|
{
|
|
t = saturate(t);
|
|
var inv = 1.0f - t;
|
|
return inv * inv * inv * p0 + 3.0f * inv * inv * t * p1 + 3.0f * inv * t * t * p2 + t * t * t * p3;
|
|
}
|
|
|
|
public static float smoothstep01(float value) => smoothstep(0.0f, 1.0f, value);
|
|
|
|
public static float smootherstep(float value)
|
|
{
|
|
value = saturate(value);
|
|
return value * value * value * (value * (value * 6.0f - 15.0f) + 10.0f);
|
|
}
|
|
|
|
public static float hash(float value) => frac(sin(value) * 43758.5453f);
|
|
public static float hash(float2 value) => hash(dot(value, new float2(127.1f, 311.7f)));
|
|
public static float hash(float3 value) => hash(dot(value, new float3(127.1f, 311.7f, 74.7f)));
|
|
|
|
public static uint asuint(float value) => (uint)BitConverter.SingleToInt32Bits(value);
|
|
public static float asfloat(uint value) => BitConverter.Int32BitsToSingle((int)value);
|
|
|
|
// PCG hashes from Jarzynski and Olano, "Hash Functions for GPU Rendering" (JCGT 9(3), 2020). Integer-only, so
|
|
// dxc and C# agree bit for bit. pcg3d/pcg4d carry uint bit patterns in int3/int4; float overloads hash IEEE bits.
|
|
public static uint pcg(uint value)
|
|
{
|
|
var state = value * 747796405u + 2891336453u;
|
|
var word = ((state >> (int)((state >> 28) + 4)) ^ state) * 277803737u;
|
|
return (word >> 22) ^ word;
|
|
}
|
|
|
|
public static int3 pcg3d(int3 value)
|
|
{
|
|
uint x = (uint)value.x * 1664525u + 1013904223u, y = (uint)value.y * 1664525u + 1013904223u, z = (uint)value.z * 1664525u + 1013904223u;
|
|
x += y * z; y += z * x; z += x * y;
|
|
x ^= x >> 16; y ^= y >> 16; z ^= z >> 16;
|
|
x += y * z; y += z * x; z += x * y;
|
|
return new int3((int)x, (int)y, (int)z);
|
|
}
|
|
|
|
// The paper hashes 2 -> N inputs through pcg3d with the unused component held constant.
|
|
public static int3 pcg3d(float2 value) => pcg3d(new float3(value, 0.0f));
|
|
public static int3 pcg3d(float3 value) => pcg3d(new int3((int)asuint(value.x), (int)asuint(value.y), (int)asuint(value.z)));
|
|
|
|
public static int4 pcg4d(int4 value)
|
|
{
|
|
uint x = (uint)value.x * 1664525u + 1013904223u, y = (uint)value.y * 1664525u + 1013904223u;
|
|
uint z = (uint)value.z * 1664525u + 1013904223u, w = (uint)value.w * 1664525u + 1013904223u;
|
|
x += y * w; y += z * x; z += x * y; w += y * z;
|
|
x ^= x >> 16; y ^= y >> 16; z ^= z >> 16; w ^= w >> 16;
|
|
x += y * w; y += z * x; z += x * y; w += y * z;
|
|
return new int4((int)x, (int)y, (int)z, (int)w);
|
|
}
|
|
|
|
public static int4 pcg4d(float4 value) =>
|
|
pcg4d(new int4((int)asuint(value.x), (int)asuint(value.y), (int)asuint(value.z), (int)asuint(value.w)));
|
|
|
|
public static float log(float value) => MathF.Log(value);
|
|
public static float2 log(float2 value) => new(log(value.x), log(value.y));
|
|
public static float3 log(float3 value) => new(log(value.x), log(value.y), log(value.z));
|
|
public static float4 log(float4 value) => new(log(value.x), log(value.y), log(value.z), log(value.w));
|
|
|
|
public static float2 pow(float2 value, float2 power) => new(pow(value.x, power.x), pow(value.y, power.y));
|
|
public static float3 pow(float3 value, float3 power) => new(pow(value.x, power.x), pow(value.y, power.y), pow(value.z, power.z));
|
|
public static float4 pow(float4 value, float4 power) => new(pow(value.x, power.x), pow(value.y, power.y), pow(value.z, power.z), pow(value.w, power.w));
|
|
|
|
public static float2 sqrt(float2 value) => new(sqrt(value.x), sqrt(value.y));
|
|
public static float3 sqrt(float3 value) => new(sqrt(value.x), sqrt(value.y), sqrt(value.z));
|
|
public static float4 sqrt(float4 value) => new(sqrt(value.x), sqrt(value.y), sqrt(value.z), sqrt(value.w));
|
|
|
|
public static float2 sin(float2 value) => new(sin(value.x), sin(value.y));
|
|
public static float3 sin(float3 value) => new(sin(value.x), sin(value.y), sin(value.z));
|
|
public static float4 sin(float4 value) => new(sin(value.x), sin(value.y), sin(value.z), sin(value.w));
|
|
|
|
public static float2 cos(float2 value) => new(cos(value.x), cos(value.y));
|
|
public static float3 cos(float3 value) => new(cos(value.x), cos(value.y), cos(value.z));
|
|
public static float4 cos(float4 value) => new(cos(value.x), cos(value.y), cos(value.z), cos(value.w));
|
|
|
|
public static float2 acos(float2 value) => new(acos(value.x), acos(value.y));
|
|
public static float3 acos(float3 value) => new(acos(value.x), acos(value.y), acos(value.z));
|
|
public static float4 acos(float4 value) => new(acos(value.x), acos(value.y), acos(value.z), acos(value.w));
|
|
|
|
public static float2 atan2(float2 y, float2 x) => new(atan2(y.x, x.x), atan2(y.y, x.y));
|
|
public static float3 atan2(float3 y, float3 x) => new(atan2(y.x, x.x), atan2(y.y, x.y), atan2(y.z, x.z));
|
|
public static float4 atan2(float4 y, float4 x) => new(atan2(y.x, x.x), atan2(y.y, x.y), atan2(y.z, x.z), atan2(y.w, x.w));
|
|
|
|
// DXIL FMin/FMax take double overloads with the same rules as float (see min/max above).
|
|
public static double min(double left, double right) => left < right || double.IsNaN(right) ? left : right;
|
|
public static double max(double left, double right) => left >= right || double.IsNaN(right) ? left : right;
|
|
|
|
public static bool any(bool2 value) => value.x || value.y;
|
|
public static bool any(bool3 value) => value.x || value.y || value.z;
|
|
public static bool any(bool4 value) => value.x || value.y || value.z || value.w;
|
|
public static bool all(bool2 value) => value.x && value.y;
|
|
public static bool all(bool3 value) => value.x && value.y && value.z;
|
|
public static bool all(bool4 value) => value.x && value.y && value.z && value.w;
|
|
|
|
// HLSL any/all on numeric vectors test components against zero (NaN != 0 counts as true).
|
|
public static bool any(float2 value) => value.x != 0.0f || value.y != 0.0f;
|
|
public static bool any(float3 value) => value.x != 0.0f || value.y != 0.0f || value.z != 0.0f;
|
|
public static bool any(float4 value) => value.x != 0.0f || value.y != 0.0f || value.z != 0.0f || value.w != 0.0f;
|
|
public static bool all(float2 value) => value.x != 0.0f && value.y != 0.0f;
|
|
public static bool all(float3 value) => value.x != 0.0f && value.y != 0.0f && value.z != 0.0f;
|
|
public static bool all(float4 value) => value.x != 0.0f && value.y != 0.0f && value.z != 0.0f && value.w != 0.0f;
|
|
public static bool any(int2 value) => value.x != 0 || value.y != 0;
|
|
public static bool any(int3 value) => value.x != 0 || value.y != 0 || value.z != 0;
|
|
public static bool any(int4 value) => value.x != 0 || value.y != 0 || value.z != 0 || value.w != 0;
|
|
public static bool all(int2 value) => value.x != 0 && value.y != 0;
|
|
public static bool all(int3 value) => value.x != 0 && value.y != 0 && value.z != 0;
|
|
public static bool all(int4 value) => value.x != 0 && value.y != 0 && value.z != 0 && value.w != 0;
|
|
|
|
// HLSL 2021 argument order: select(condition, whenTrue, whenFalse).
|
|
// Unity.Mathematics' select(falseValue, trueValue, condition) is the reverse.
|
|
public static float select(bool condition, float whenTrue, float whenFalse) => condition ? whenTrue : whenFalse;
|
|
public static float2 select(bool2 condition, float2 whenTrue, float2 whenFalse) =>
|
|
new(condition.x ? whenTrue.x : whenFalse.x, condition.y ? whenTrue.y : whenFalse.y);
|
|
public static float3 select(bool3 condition, float3 whenTrue, float3 whenFalse) =>
|
|
new(condition.x ? whenTrue.x : whenFalse.x, condition.y ? whenTrue.y : whenFalse.y, condition.z ? whenTrue.z : whenFalse.z);
|
|
public static float4 select(bool4 condition, float4 whenTrue, float4 whenFalse) =>
|
|
new(condition.x ? whenTrue.x : whenFalse.x, condition.y ? whenTrue.y : whenFalse.y, condition.z ? whenTrue.z : whenFalse.z, condition.w ? whenTrue.w : whenFalse.w);
|
|
|
|
// HLSL mul: mul(m, v) treats v as a column vector, mul(v, m) as a row vector.
|
|
public static float mul(float left, float right) => left * right;
|
|
|
|
public static float2 mul(float2x2 m, float2 v) =>
|
|
new(m._m00 * v.x + m._m01 * v.y, m._m10 * v.x + m._m11 * v.y);
|
|
|
|
public static float2 mul(float2 v, float2x2 m) =>
|
|
new(v.x * m._m00 + v.y * m._m10, v.x * m._m01 + v.y * m._m11);
|
|
|
|
public static float2x2 mul(float2x2 a, float2x2 b) =>
|
|
new(
|
|
a._m00 * b._m00 + a._m01 * b._m10, a._m00 * b._m01 + a._m01 * b._m11,
|
|
a._m10 * b._m00 + a._m11 * b._m10, a._m10 * b._m01 + a._m11 * b._m11);
|
|
|
|
public static float3 mul(float3x3 m, float3 v) =>
|
|
new(
|
|
m._m00 * v.x + m._m01 * v.y + m._m02 * v.z,
|
|
m._m10 * v.x + m._m11 * v.y + m._m12 * v.z,
|
|
m._m20 * v.x + m._m21 * v.y + m._m22 * v.z);
|
|
|
|
public static float3 mul(float3 v, float3x3 m) =>
|
|
new(
|
|
v.x * m._m00 + v.y * m._m10 + v.z * m._m20,
|
|
v.x * m._m01 + v.y * m._m11 + v.z * m._m21,
|
|
v.x * m._m02 + v.y * m._m12 + v.z * m._m22);
|
|
|
|
public static float3x3 mul(float3x3 a, float3x3 b) =>
|
|
new(
|
|
a._m00 * b._m00 + a._m01 * b._m10 + a._m02 * b._m20,
|
|
a._m00 * b._m01 + a._m01 * b._m11 + a._m02 * b._m21,
|
|
a._m00 * b._m02 + a._m01 * b._m12 + a._m02 * b._m22,
|
|
a._m10 * b._m00 + a._m11 * b._m10 + a._m12 * b._m20,
|
|
a._m10 * b._m01 + a._m11 * b._m11 + a._m12 * b._m21,
|
|
a._m10 * b._m02 + a._m11 * b._m12 + a._m12 * b._m22,
|
|
a._m20 * b._m00 + a._m21 * b._m10 + a._m22 * b._m20,
|
|
a._m20 * b._m01 + a._m21 * b._m11 + a._m22 * b._m21,
|
|
a._m20 * b._m02 + a._m21 * b._m12 + a._m22 * b._m22);
|
|
|
|
// Abramowitz & Stegun 7.1.26, a single-precision rational/exponential fit with a stated
|
|
// max absolute error of 1.5e-7. erf is odd, so only x >= 0 is fit and the sign is restored.
|
|
// Measured against an independent reference: worst absolute error 6.621e-7, inside the
|
|
// 2e-6 bar this package's tests pin.
|
|
public static float erf(float value)
|
|
{
|
|
const float p = 0.3275911f;
|
|
const float a1 = 0.254829592f;
|
|
const float a2 = -0.284496736f;
|
|
const float a3 = 1.421413741f;
|
|
const float a4 = -1.453152027f;
|
|
const float a5 = 1.061405429f;
|
|
|
|
var s = sign(value);
|
|
var x = abs(value);
|
|
var t = 1.0f / (1.0f + p * x);
|
|
var poly = ((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t;
|
|
var result = 1.0f - poly * exp(-x * x);
|
|
return s * result;
|
|
}
|
|
|
|
// Giles, "Approximating the erfinv function" (GPU Computing Gems, 2010): a single-precision
|
|
// minimax polynomial in two ranges of w = -ln((1-x)(1+x)), central and tail.
|
|
// Domain is [-1, 1]. Outside it, including +-infinity, the result is NaN (there is no real
|
|
// inverse). At the edges the mathematical limit is infinite: erfinv(1) = +infinity and
|
|
// erfinv(-1) = -infinity. The unguarded polynomial hits log(0) = -infinity at those edges and
|
|
// the tail branch's leading coefficient is negative, so it delivers the wrong sign; guard the
|
|
// edges explicitly instead of trusting the fit through the singularity.
|
|
// Measured against an independent reference: worst absolute error 5.066e-7, and the
|
|
// erf(erfinv(y)) round trip 6.109e-7 worst case, both inside the 2e-6 bar this package's
|
|
// tests pin.
|
|
public static float erfinv(float value)
|
|
{
|
|
if (value >= 1.0f) return value > 1.0f ? float.NaN : float.PositiveInfinity;
|
|
if (value <= -1.0f) return value < -1.0f ? float.NaN : float.NegativeInfinity;
|
|
|
|
var w = -log((1.0f - value) * (1.0f + value));
|
|
float p;
|
|
if (w < 5.0f)
|
|
{
|
|
w -= 2.5f;
|
|
p = 2.81022636e-08f;
|
|
p = 3.43273939e-07f + p * w;
|
|
p = -3.5233877e-06f + p * w;
|
|
p = -4.39150654e-06f + p * w;
|
|
p = 0.00021858087f + p * w;
|
|
p = -0.00125372503f + p * w;
|
|
p = -0.00417768164f + p * w;
|
|
p = 0.246640727f + p * w;
|
|
p = 1.50140941f + p * w;
|
|
}
|
|
else
|
|
{
|
|
w = sqrt(w) - 3.0f;
|
|
p = -0.000200214257f;
|
|
p = 0.000100950558f + p * w;
|
|
p = 0.00134934322f + p * w;
|
|
p = -0.00367342844f + p * w;
|
|
p = 0.00573950773f + p * w;
|
|
p = -0.0076224613f + p * w;
|
|
p = 0.00943887047f + p * w;
|
|
p = 1.00167406f + p * w;
|
|
p = 2.83297682f + p * w;
|
|
}
|
|
|
|
return p * value;
|
|
}
|
|
|
|
// Ashima Arts / Ian McEwan 3D simplex noise (MIT), mirrored by cultmath_snoise(float3)
|
|
// in shaders/CultMath.hlsl with the same float32 evaluation order. Follows stegu/webgl-noise
|
|
// src/noise3D.glsl at 22434e04d7 (kernel radius^2 0.5, scale 105, permute (34x+10)x); see
|
|
// THIRD-PARTY-NOTICES.md. Continuous across simplex-cell boundaries
|
|
// (NoiseGradTests.SnoiseAndItsGradientAreContinuousAcrossSimplexCellBoundaries). Ties in x0 resolve
|
|
// in Gustavson's simplexnoise1234 order (x over y, y over z, x over z), a total order, so the four
|
|
// corners are always a simplex; upstream's step is cyclic and degenerates where x0's components are
|
|
// equal (NoiseGradTests.SnoiseIsContinuousOnTheSimplexDiagonal).
|
|
// SNOISE_LIPSCHITZ (math.Interval.cs) is measured from this kernel: a change here must re-run
|
|
// NoiseBoundTests.MeasureLipschitz and re-pin it.
|
|
public static float snoise(float3 value)
|
|
{
|
|
const float cx = 1.0f / 6.0f;
|
|
const float cy = 1.0f / 3.0f;
|
|
var i = floor(value + dot(value, new float3(cy, cy, cy)));
|
|
var x0 = value - i + dot(i, new float3(cx, cx, cx));
|
|
var g = new float3(step(x0.y, x0.x), step(x0.z, x0.y), 1.0f - step(x0.z, x0.x));
|
|
var l = 1.0f - g;
|
|
var lzxy = new float3(l.z, l.x, l.y);
|
|
var i1 = min(g, lzxy);
|
|
var i2 = max(g, lzxy);
|
|
var x1 = x0 - i1 + cx;
|
|
var x2 = x0 - i2 + cy;
|
|
var x3 = x0 - 0.5f;
|
|
|
|
i = snoise_mod289(i);
|
|
var p = snoise_permute(snoise_permute(snoise_permute(
|
|
i.z + new float4(0.0f, i1.z, i2.z, 1.0f)) + i.y + new float4(0.0f, i1.y, i2.y, 1.0f)) + i.x + new float4(0.0f, i1.x, i2.x, 1.0f));
|
|
|
|
const float n = 0.142857142857f;
|
|
var ns = new float3(2.0f * n, 0.5f * n - 1.0f, n);
|
|
var j = p - 49.0f * floor(p * ns.z * ns.z);
|
|
var xs = floor(j * ns.z);
|
|
var ys = floor(j - 7.0f * xs);
|
|
var x = xs * ns.x + ns.y;
|
|
var y = ys * ns.x + ns.y;
|
|
var h = 1.0f - abs(x) - abs(y);
|
|
|
|
var b0 = new float4(x.x, x.y, y.x, y.y);
|
|
var b1 = new float4(x.z, x.w, y.z, y.w);
|
|
var s0 = floor(b0) * 2.0f + 1.0f;
|
|
var s1 = floor(b1) * 2.0f + 1.0f;
|
|
var sh = -step(h, 0.0f);
|
|
var a0 = new float4(b0.x, b0.z, b0.y, b0.w) + new float4(s0.x, s0.z, s0.y, s0.w) * new float4(sh.x, sh.x, sh.y, sh.y);
|
|
var a1 = new float4(b1.x, b1.z, b1.y, b1.w) + new float4(s1.x, s1.z, s1.y, s1.w) * new float4(sh.z, sh.z, sh.w, sh.w);
|
|
|
|
var p0 = new float3(a0.x, a0.y, h.x);
|
|
var p1 = new float3(a0.z, a0.w, h.y);
|
|
var p2 = new float3(a1.x, a1.y, h.z);
|
|
var p3 = new float3(a1.z, a1.w, h.w);
|
|
var norm = 1.79284291400159f - 0.85373472095314f * new float4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3));
|
|
p0 *= norm.x;
|
|
p1 *= norm.y;
|
|
p2 *= norm.z;
|
|
p3 *= norm.w;
|
|
|
|
var m = max(0.5f - new float4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), 0.0f);
|
|
m = m * m;
|
|
return 105.0f * dot(m * m, new float4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3)));
|
|
}
|
|
|
|
// Ashima Arts / Ian McEwan 2D simplex noise, kept component-explicit so
|
|
// the C# and HLSL mirrors preserve the same float32 evaluation order.
|
|
public static float snoise(float2 value)
|
|
{
|
|
var c = new float4(
|
|
0.211324865405187f,
|
|
0.366025403784439f,
|
|
-0.577350269189626f,
|
|
0.024390243902439f);
|
|
var i = floor(value + dot(value, new float2(c.y, c.y)));
|
|
var x0 = value - i + dot(i, new float2(c.x, c.x));
|
|
var i1 = x0.x > x0.y ? new float2(1.0f, 0.0f) : new float2(0.0f, 1.0f);
|
|
var x12 = new float4(
|
|
x0.x + c.x - i1.x,
|
|
x0.y + c.x - i1.y,
|
|
x0.x + c.z,
|
|
x0.y + c.z);
|
|
|
|
i = snoise_mod289(i);
|
|
var p = snoise_permute(snoise_permute(i.y + new float3(0.0f, i1.y, 1.0f)) + i.x + new float3(0.0f, i1.x, 1.0f));
|
|
var m = max(0.5f - new float3(
|
|
dot(x0, x0),
|
|
x12.x * x12.x + x12.y * x12.y,
|
|
x12.z * x12.z + x12.w * x12.w), 0.0f);
|
|
m = m * m;
|
|
m = m * m;
|
|
|
|
var x = 2.0f * frac(p * c.w) - 1.0f;
|
|
var h = abs(x) - 0.5f;
|
|
var ox = floor(x + 0.5f);
|
|
var a0 = x - ox;
|
|
m *= 1.79284291400159f - 0.85373472095314f * (a0 * a0 + h * h);
|
|
|
|
var g = new float3(
|
|
a0.x * x0.x + h.x * x0.y,
|
|
a0.y * x12.x + h.y * x12.y,
|
|
a0.z * x12.z + h.z * x12.w);
|
|
return 130.0f * dot(m, g);
|
|
}
|
|
|
|
private static float2 snoise_mod289(float2 value) => value - floor(value * (1.0f / 289.0f)) * 289.0f;
|
|
private static float3 snoise_mod289(float3 value) => value - floor(value * (1.0f / 289.0f)) * 289.0f;
|
|
private static float4 snoise_mod289(float4 value) => value - floor(value * (1.0f / 289.0f)) * 289.0f;
|
|
private static float3 snoise_permute(float3 value) => snoise_mod289(((value * 34.0f) + 10.0f) * value);
|
|
private static float4 snoise_permute(float4 value) => snoise_mod289(((value * 34.0f) + 10.0f) * value);
|
|
|
|
// Analytic gradient of snoise(float3) (invariant 8), following the differentiation in
|
|
// webgl-noise src/noise3Dgrad.glsl (Ashima Arts / McEwan, MIT), which uses the same kernel
|
|
// radius^2 0.5, permutation (34x+10)x and scale 105 as snoise(float3) (stegu/webgl-noise
|
|
// 22434e04d7). Every local up through p0..p3 and x0..x3 below is snoise(float3)'s own derivation
|
|
// verbatim, so .w is that same field (NoiseGradTests.SnoiseGradValueIsBitEqualToSnoise pins the
|
|
// two bit-equal over a random sweep). Per corner, m0 = max(0.5 - dot(x,x), 0), so dm0/dx = -2x,
|
|
// and d(m0^4 * dot(p,x))/dx = 4*m0^3*(-2x)*dot(p,x) + m0^4*p = -8*m0^3*dot(p,x)*x + m0^4*p (p is
|
|
// that corner's own gradient constant, independent of x); summed over the four corners and
|
|
// scaled by the same 105 as the value.
|
|
public static float4 snoise_grad(float3 value)
|
|
{
|
|
const float cx = 1.0f / 6.0f;
|
|
const float cy = 1.0f / 3.0f;
|
|
var i = floor(value + dot(value, new float3(cy, cy, cy)));
|
|
var x0 = value - i + dot(i, new float3(cx, cx, cx));
|
|
var g = new float3(step(x0.y, x0.x), step(x0.z, x0.y), 1.0f - step(x0.z, x0.x));
|
|
var l = 1.0f - g;
|
|
var lzxy = new float3(l.z, l.x, l.y);
|
|
var i1 = min(g, lzxy);
|
|
var i2 = max(g, lzxy);
|
|
var x1 = x0 - i1 + cx;
|
|
var x2 = x0 - i2 + cy;
|
|
var x3 = x0 - 0.5f;
|
|
|
|
i = snoise_mod289(i);
|
|
var p = snoise_permute(snoise_permute(snoise_permute(
|
|
i.z + new float4(0.0f, i1.z, i2.z, 1.0f)) + i.y + new float4(0.0f, i1.y, i2.y, 1.0f)) + i.x + new float4(0.0f, i1.x, i2.x, 1.0f));
|
|
|
|
const float n = 0.142857142857f;
|
|
var ns = new float3(2.0f * n, 0.5f * n - 1.0f, n);
|
|
var j = p - 49.0f * floor(p * ns.z * ns.z);
|
|
var xs = floor(j * ns.z);
|
|
var ys = floor(j - 7.0f * xs);
|
|
var x = xs * ns.x + ns.y;
|
|
var y = ys * ns.x + ns.y;
|
|
var h = 1.0f - abs(x) - abs(y);
|
|
|
|
var b0 = new float4(x.x, x.y, y.x, y.y);
|
|
var b1 = new float4(x.z, x.w, y.z, y.w);
|
|
var s0 = floor(b0) * 2.0f + 1.0f;
|
|
var s1 = floor(b1) * 2.0f + 1.0f;
|
|
var sh = -step(h, 0.0f);
|
|
var a0 = new float4(b0.x, b0.z, b0.y, b0.w) + new float4(s0.x, s0.z, s0.y, s0.w) * new float4(sh.x, sh.x, sh.y, sh.y);
|
|
var a1 = new float4(b1.x, b1.z, b1.y, b1.w) + new float4(s1.x, s1.z, s1.y, s1.w) * new float4(sh.z, sh.z, sh.w, sh.w);
|
|
|
|
var p0 = new float3(a0.x, a0.y, h.x);
|
|
var p1 = new float3(a0.z, a0.w, h.y);
|
|
var p2 = new float3(a1.x, a1.y, h.z);
|
|
var p3 = new float3(a1.z, a1.w, h.w);
|
|
var norm = 1.79284291400159f - 0.85373472095314f * new float4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3));
|
|
p0 *= norm.x;
|
|
p1 *= norm.y;
|
|
p2 *= norm.z;
|
|
p3 *= norm.w;
|
|
|
|
var m0 = max(0.5f - new float4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), 0.0f);
|
|
var m2 = m0 * m0;
|
|
var m3 = m2 * m0;
|
|
var m4 = m2 * m2;
|
|
|
|
var px = new float4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3));
|
|
var value2 = 105.0f * dot(m4, px);
|
|
var grad = -8.0f * m3.x * x0 * px.x + m4.x * p0
|
|
+ -8.0f * m3.y * x1 * px.y + m4.y * p1
|
|
+ -8.0f * m3.z * x2 * px.z + m4.z * p2
|
|
+ -8.0f * m3.w * x3 * px.w + m4.w * p3;
|
|
grad *= 105.0f;
|
|
|
|
return new float4(grad, value2);
|
|
}
|
|
|
|
// Octave sum of snoise_grad (invariant 8). Each octave samples snoise_grad at p scaled by the
|
|
// running frequency, so by the chain rule its gradient scales by that same frequency
|
|
// (d/dp[n(f*p)] = f * grad_n(f*p)); amplitude and frequency both compound per octave (gain,
|
|
// lacunarity respectively), not applied once for the whole sum, so with 1 octave this is exactly
|
|
// snoise_grad (NoiseGradTests.FbmGradWithOneOctaveEqualsSnoiseGrad). octaves is clamped to
|
|
// [0, 16]: HlslSourceCompatibilityTests' bit-parity oracle drives every mirrored int parameter
|
|
// across the full int32 range, which is the right domain for pcg3d/pcg4d's O(1) hash inputs but
|
|
// would turn an unclamped octave count into a multi-billion-iteration loop for this mirror.
|
|
public static float4 fbm_grad(float3 p, int octaves, float lacunarity, float gain)
|
|
{
|
|
octaves = clamp(octaves, 0, 16);
|
|
var amplitude = 1.0f;
|
|
var frequency = 1.0f;
|
|
var value = 0.0f;
|
|
var gradient = new float3(0.0f, 0.0f, 0.0f);
|
|
for (var i = 0; i < octaves; i++)
|
|
{
|
|
var n = snoise_grad(p * frequency);
|
|
value += amplitude * n.w;
|
|
gradient += amplitude * frequency * new float3(n.x, n.y, n.z);
|
|
frequency *= lacunarity;
|
|
amplitude *= gain;
|
|
}
|
|
|
|
return new float4(gradient, value);
|
|
}
|
|
|
|
// Ridged multifractal noise (Musgrave, Texturing and Modeling, ch. 16; the operator's dune term,
|
|
// 2026-09-25): each octave folds snoise_grad about zero, value = Sum a_i*(1 - |n_i|), gradient =
|
|
// -Sum a_i*f_i*sign(n_i)*grad(n_i), the same per-octave frequency scaling as fbm_grad's chain
|
|
// rule plus the fold's sign flip (d|n|/dp = sign(n)*dn/dp). The crease at n_i = 0 is real and
|
|
// deliberate, not a bug: NoiseGradTests excludes a band around it, the same way cellular excludes
|
|
// its F1 = F2 seam. sign follows this codebase's HLSL-matching convention (design.md: sign
|
|
// returns int, 0 at exactly 0), so the crease itself contributes no gradient rather than an
|
|
// arbitrary one. Same octave clamp as fbm_grad, for the same bit-parity-oracle reason.
|
|
public static float4 ridged_grad(float3 p, int octaves, float lacunarity, float gain)
|
|
{
|
|
octaves = clamp(octaves, 0, 16);
|
|
var amplitude = 1.0f;
|
|
var frequency = 1.0f;
|
|
var value = 0.0f;
|
|
var gradient = new float3(0.0f, 0.0f, 0.0f);
|
|
for (var i = 0; i < octaves; i++)
|
|
{
|
|
var n = snoise_grad(p * frequency);
|
|
var s = (float)sign(n.w);
|
|
value += amplitude * (1.0f - abs(n.w));
|
|
gradient += -amplitude * frequency * s * new float3(n.x, n.y, n.z);
|
|
frequency *= lacunarity;
|
|
amplitude *= gain;
|
|
}
|
|
|
|
return new float4(gradient, value);
|
|
}
|
|
|
|
// Inigo Quilez, "smooth minimum" (https://iquilezles.org/articles/smin/): the quadratic-polynomial
|
|
// smin, carried to value-and-gradient form (invariant 8). h is affine in (b.w - a.w) inside the
|
|
// smoothing band, so differentiating value = lerp(b.w, a.w, h) - k*h*(1-h) with respect to
|
|
// position, the terms carrying dh/dp cancel exactly (dh/dp * [(a.w-b.w) - k + 2*k*h] and the
|
|
// bracket is identically zero given how h is built from a.w-b.w). What is left is exactly
|
|
// lerp(∇b, ∇a, h): the gradient blend is not an approximation, it is the analytic gradient.
|
|
// Outside the band (|a.w - b.w| >= k), h saturates to 0 or 1 and this is exactly min(a, b) with
|
|
// that input's own gradient, continuously (h and the correction term both reach the boundary at
|
|
// the same value from the smooth side, so there is no kink to exclude here the way cellular's
|
|
// F1 = F2 seam needs one).
|
|
public static float4 smin_grad(float4 a, float4 b, float k)
|
|
{
|
|
var h = saturate(0.5f + 0.5f * (b.w - a.w) / k);
|
|
var value = lerp(b.w, a.w, h) - k * h * (1.0f - h);
|
|
var gradient = lerp(new float3(b.x, b.y, b.z), new float3(a.x, a.y, a.z), h);
|
|
return new float4(gradient, value);
|
|
}
|
|
|
|
// Worley, "A Cellular Texture Basis Function" (SIGGRAPH 1996): F1/F2 and their analytic
|
|
// gradients, searched over the jittered neighbourhood of feature points. Feature points are
|
|
// selected with pcg3d over each neighbour's integer cell coordinate (design.md, "Integer
|
|
// hashing"), never the float32 bit pattern of that coordinate and never the sin-based hash.
|
|
// Hashing the integer cell coordinate directly is exact and runtime-independent, and it is what
|
|
// the HLSL mirror does, so C# and shader stay bit-for-bit aligned. What pins this is the oracle
|
|
// and mirror tests, not a uniformity measurement: reverting to hashing the float32 bit pattern
|
|
// (or the sin-based hash) is caught by ProductionSearchMatchesTheBruteForceOracleExactly,
|
|
// EveryMirrorFunctionMatchesCSharpMath, F1NeverExceedsF2, and the other tests that key off the
|
|
// integer-hash fixtures. The
|
|
// gradient of a distance field, del|p - c|, is the unit vector (p - c)/|p - c|; that is undefined exactly at
|
|
// a feature point (F1 = 0 or F2 = 0), so cellular follows the repo's existing degenerate-normal
|
|
// convention (GameCult.Geometry.CultGeometryIsoSurface.EmitOrientedTriangle: guard the zero-length
|
|
// case and return the zero vector instead of the NaN a bare normalize would produce there).
|
|
//
|
|
// The search radius is 2 cells on every axis (5x5x5), not the classic 1-cell (3x3x3) radius,
|
|
// because 3x3x3 is not exact: measured 9 of 400,000 points in [-50,50]^3 land on the wrong F2
|
|
// (worst error 0.088), and past |p| of about 2^23 the float32 jitter can put two distinct cells'
|
|
// features on the same value, which 3x3x3 has no chance of noticing either way. 5x5x5 is exact,
|
|
// not merely safer, by the following argument:
|
|
//
|
|
// Let u = p - cell be the query's position inside its own cell (u in [0,1)^3 per component), and
|
|
// let a candidate feature at integer cell offset d = (dx, dy, dz) sit at d + j for some jitter j
|
|
// in [0,1)^3. Per axis, the displacement from p to that feature is (u_i - j_i) - d_i; since
|
|
// u_i - j_i is always strictly inside (-1, 1) regardless of the actual u_i and j_i values,
|
|
// |displacement_i| is always strictly greater than max(0, |d_i| - 1). So every cell with |d_i| >=
|
|
// 3 on some axis is strictly farther than 2 from p (that axis alone already gives max(0, 3-1) =
|
|
// 2). Now bound F1 and F2 from above using two real, always-present candidates: the query's own
|
|
// cell (d = 0), whose distance is always strictly less than sqrt(3) (every axis strictly under
|
|
// 1); and the near-side neighbour on whichever single axis's u_i sits closest to a cell
|
|
// boundary. Let m = min(u_i, 1 - u_i) on that chosen axis (m is in [0, 1/2] by construction,
|
|
// since it is the smaller of the two). Stepping one cell towards that nearer boundary puts the
|
|
// query within (1 + m) of the far face on the offset axis and within (1 - m) of the near face on
|
|
// each of the other two axes, so that neighbour's squared distance is at most
|
|
// (1 + m)^2 + 2*(1 - m)^2 = 3 - 2m + 3m^2, which is <= 3 for every m in [0, 1/2] (3 at m = 0,
|
|
// dipping to a minimum of 8/3 ~= 2.667 at m = 1/3, then rising back to 2.75 at m = 1/2), and
|
|
// strictly below 3 because the jitter never reaches exactly 0 or
|
|
// 1 (u_i - j_i stays strictly inside (-1, 1), so this bound, like the own-cell one, is never
|
|
// tight). F2 is at most the larger of these two real candidates (adding candidates to a set can
|
|
// only lower or hold its 2nd-smallest value), so F1 < sqrt(3) and F2 < sqrt(3) always: roughly
|
|
// 1.732 < 2, margin 2 - sqrt(3) ~= 0.268 (Soul's adversarial check, cut 2a-i: the worst
|
|
// realizable in-box F2 is exactly sqrt(3), reached only in the unattainable limit at a cell
|
|
// corner). Every cell outside the 5x5x5 block is farther than 2, hence farther than either F1 or
|
|
// F2, so the true global F1 and F2 are always inside it. (The matching argument for 3x3x3 only
|
|
// reaches radius 1, whose excluded cells start at distance > 1, which is less than the 1.732
|
|
// ceiling — the gap 3x3x3's failures live in.)
|
|
//
|
|
// An earlier version of this proof picked the near-side neighbour on a fixed axis rather than
|
|
// whichever axis sits closest to its boundary, and that version is wrong (Soul, cut 2a-i):
|
|
// u = (0.5, 0, 0) with the x-neighbour's jittered feature at (-1, 1 - 2^-24, 1 - 2^-24) gives a
|
|
// distance of about 2.06, outside the radius-2 exclusion the proof needs. Choosing the axis by
|
|
// proximity to its own boundary is what keeps m <= 1/2 and the bound at sqrt(3); x was the wrong
|
|
// axis to fix there because u.x = 0.5 is the point on the whole cube farthest from any boundary.
|
|
//
|
|
// The loop below still visits every one of the 125 cells in a fixed order, but skips the
|
|
// pcg3d/jitter/distance work for a cell whose distance cannot possibly beat the current f2: that
|
|
// per-cell lower bound (the max(0, |d_i| - 1) term above, combined across axes) depends only on
|
|
// the integer offset, never on p or the hash, so pruning a cell the unpruned loop would also
|
|
// have rejected changes no output bit (CellularExactSearchTests pins the pruned/unpruned
|
|
// equivalence and the exactness against a 7x7x7 brute-force oracle).
|
|
public static CultCellular cellular(float3 p)
|
|
{
|
|
var cell = floor(p);
|
|
var f1 = 1.0e30f;
|
|
var f2 = 1.0e30f;
|
|
var c1 = new float3(0.0f, 0.0f, 0.0f);
|
|
var c2 = new float3(0.0f, 0.0f, 0.0f);
|
|
var cellId = new float3(0.0f, 0.0f, 0.0f);
|
|
|
|
for (var dz = -2; dz <= 2; dz++)
|
|
{
|
|
var lz = max(0.0f, abs((float)dz) - 1.0f);
|
|
for (var dy = -2; dy <= 2; dy++)
|
|
{
|
|
var ly = max(0.0f, abs((float)dy) - 1.0f);
|
|
for (var dx = -2; dx <= 2; dx++)
|
|
{
|
|
var lx = max(0.0f, abs((float)dx) - 1.0f);
|
|
var lowerBound = sqrt(lx * lx + ly * ly + lz * lz);
|
|
if (lowerBound >= f2)
|
|
continue;
|
|
|
|
var neighbor = cell + new float3(dx, dy, dz);
|
|
var hash = pcg3d(int3(neighbor));
|
|
var jitter = new float3(cellular_unit(hash.x), cellular_unit(hash.y), cellular_unit(hash.z));
|
|
var feature = neighbor + jitter;
|
|
var d = distance(p, feature);
|
|
|
|
if (d < f1)
|
|
{
|
|
f2 = f1; c2 = c1;
|
|
f1 = d; c1 = feature; cellId = neighbor;
|
|
}
|
|
else if (d < f2)
|
|
{
|
|
f2 = d; c2 = feature;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// F1 = 0 is reachable (p sits exactly on a feature point). F2 = 0 needs two distinct cells'
|
|
// pcg3d-jittered features to land on the same float32 value in every component: unreachable
|
|
// at ordinary magnitudes, but confirmed reachable past roughly |p| = 2^23 (design.md,
|
|
// "Precision domain"), so both share the same degenerate-gradient guard.
|
|
var grad1 = f1 > 0.0f ? (p - c1) / f1 : new float3(0.0f, 0.0f, 0.0f);
|
|
var grad2 = f2 > 0.0f ? (p - c2) / f2 : new float3(0.0f, 0.0f, 0.0f);
|
|
|
|
// id comes from pcg4d over the winning cell's integer coordinate alone (same int-hashing
|
|
// rule as the jitter above), never from the pcg3d hash that produced its jitter: sharing a
|
|
// hash between id and jitter would make id a deterministic function of the jitter it is
|
|
// supposed to be independent from, correlated across every cell (F4).
|
|
var idHash = pcg4d(new int4(int3(cellId), 0));
|
|
|
|
return new CultCellular(
|
|
new float4(grad1, f1),
|
|
new float4(grad2 - grad1, f2 - f1),
|
|
cellular_unit(idHash.w));
|
|
}
|
|
|
|
// Maps a hash component's uint bit pattern to [0, 1) using the top 24 bits, so the result is
|
|
// always strictly less than 1: a plain divide-by-2^32 rounds 0xFFFFFF80..0xFFFFFFFF up to
|
|
// exactly 1.0f in float32, outside the documented range. (uint >> 8) has at most 24 significant
|
|
// bits, which float32's 24-bit mantissa represents exactly, and the division is by an exact power
|
|
// of two, so both the shift-and-convert and the divide are the same IEEE-754 round-to-nearest on
|
|
// dxc and C#, bit-exact on both sides.
|
|
internal static float cellular_unit(int bits) => ((uint)bits >> 8) * (1.0f / 16777216.0f);
|
|
|
|
public static float value_noise(float2 position)
|
|
{
|
|
var cell = floor(position);
|
|
var local = frac(position);
|
|
var u = local * local * (3.0f - 2.0f * local);
|
|
|
|
var a = hash(cell);
|
|
var b = hash(cell + new float2(1.0f, 0.0f));
|
|
var c = hash(cell + new float2(0.0f, 1.0f));
|
|
var d = hash(cell + new float2(1.0f, 1.0f));
|
|
|
|
return lerp(lerp(a, b, u.x), lerp(c, d, u.x), u.y);
|
|
}
|
|
|
|
public static float value_noise_bicubic(float2 position)
|
|
{
|
|
var cell = floor(position);
|
|
var local = frac(position);
|
|
|
|
var y0 = catmullrom(
|
|
hash(cell + new float2(-1.0f, -1.0f)),
|
|
hash(cell + new float2(0.0f, -1.0f)),
|
|
hash(cell + new float2(1.0f, -1.0f)),
|
|
hash(cell + new float2(2.0f, -1.0f)),
|
|
local.x);
|
|
var y1 = catmullrom(
|
|
hash(cell + new float2(-1.0f, 0.0f)),
|
|
hash(cell + new float2(0.0f, 0.0f)),
|
|
hash(cell + new float2(1.0f, 0.0f)),
|
|
hash(cell + new float2(2.0f, 0.0f)),
|
|
local.x);
|
|
var y2 = catmullrom(
|
|
hash(cell + new float2(-1.0f, 1.0f)),
|
|
hash(cell + new float2(0.0f, 1.0f)),
|
|
hash(cell + new float2(1.0f, 1.0f)),
|
|
hash(cell + new float2(2.0f, 1.0f)),
|
|
local.x);
|
|
var y3 = catmullrom(
|
|
hash(cell + new float2(-1.0f, 2.0f)),
|
|
hash(cell + new float2(0.0f, 2.0f)),
|
|
hash(cell + new float2(1.0f, 2.0f)),
|
|
hash(cell + new float2(2.0f, 2.0f)),
|
|
local.x);
|
|
|
|
return catmullrom(y0, y1, y2, y3, local.y);
|
|
}
|
|
}
|