From 7904e15ed39b72b0b5ce79443019e51bdc99ad44 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 23:51:02 +0200 Subject: [PATCH 1/4] CultMath: phacelle, Johansen's Phacelle noise on 3D cells with an exact gradient Cut 2a-ii-b. phacelle(float3 p, float3 side, float offset, float normalization) returns CultPhasor { float4 cos; float4 sin }, each float4(gradient.xyz, value.w). - 4x4x4 cells, pcg3d-hashed jitter in [-0.5, 0.5), weight max(0, exp(-2 d^2) - 0.01111), normalization as upstream (lpmitchell/AdvancedTerrainErosion PhacelleNoise). The caller supplies the stripe wave vector; offset is in cycles. - The gradient is exact: weight derivatives and the normalization are included, not upstream's grad(cos) ~ -sin * side. - Cells whose per-axis lower bound is at least 2.25 have weight exactly 0 and are skipped. An internal overload switches the prune off so the test can pin bit equality without a second copy of the loop. - cultmath_phacelle and CultPhasor mirror it in both tracked CultMath.hlsl copies (the same blob). The mirror test compares the struct field by field; it now covers 34 functions. - PhacelleTests: value against an independent double-precision reference, quarter-cycle phase rotation, gradient against central differences on both sides of the normalization threshold (excluding the threshold and each cell's weight-support edge, where the exact gradient is discontinuous), stripes along the flow on the unit sphere (0.215 against the 0.25 bound), step-difference continuity across cell planes, and pruned == unpruned bit for bit. Co-Authored-By: Claude Opus 5.5 --- packages/cultmath/shaders/CultMath.hlsl | 86 ++++++ packages/cultmath/src/CultMath/CultPhasor.cs | 21 ++ .../cultmath/src/CultMath/CultPhasor.cs.meta | 2 + packages/cultmath/src/CultMath/math.cs | 98 +++++++ .../HlslSourceCompatibilityTests.cs | 3 +- .../tests/CultMath.Tests/PhacelleTests.cs | 275 ++++++++++++++++++ .../Shaders/CultMath.hlsl | 86 ++++++ 7 files changed, 570 insertions(+), 1 deletion(-) create mode 100644 packages/cultmath/src/CultMath/CultPhasor.cs create mode 100644 packages/cultmath/src/CultMath/CultPhasor.cs.meta create mode 100644 packages/cultmath/tests/CultMath.Tests/PhacelleTests.cs diff --git a/packages/cultmath/shaders/CultMath.hlsl b/packages/cultmath/shaders/CultMath.hlsl index 50c5a746..d0dd22ae 100644 --- a/packages/cultmath/shaders/CultMath.hlsl +++ b/packages/cultmath/shaders/CultMath.hlsl @@ -503,6 +503,92 @@ CultCellular cultmath_cellular(float3 p) return result; } +// Invariant 8's struct return shape for cultmath_phacelle, compared field by field and bit for bit +// against CultMath.CultPhasor, same field order. +struct CultPhasor +{ + float4 cos; + float4 sin; +}; + +// Johansen's Phacelle noise on 3D cells. See math.phacelle's comment (math.cs) for the stripe wave +// vector, the weight, the normalization, the exact gradient and why the prune is exact. Mirrors +// the C# body operation for operation; the C#-only internal overload that switches the prune off +// has no shader twin. +CultPhasor cultmath_phacelle(float3 p, float3 side, float offset, float normalization) +{ + float3 cell = floor(p); + float3 local = p - cell; + float phaseOffset = offset * CULTMATH_TAU; + float sumW = 0.0; + float sumCos = 0.0; + float sumSin = 0.0; + float3 gradW = float3(0.0, 0.0, 0.0); + float3 gradCos = float3(0.0, 0.0, 0.0); + float3 gradSin = float3(0.0, 0.0, 0.0); + + for (int dz = -1; dz <= 2; dz++) + { + float gz = max(0.0, abs(local.z - (float)dz) - 0.5); + for (int dy = -1; dy <= 2; dy++) + { + float gy = max(0.0, abs(local.y - (float)dy) - 0.5); + for (int dx = -1; dx <= 2; dx++) + { + float gx = max(0.0, abs(local.x - (float)dx) - 0.5); + if (gx * gx + gy * gy + gz * gz >= 2.25) + continue; + + float3 gridStep = float3(dx, dy, dz); + int3 hash = cultmath_pcg3d(int3(cell + gridStep)); + float3 jitter = float3( + ((uint)hash.x >> 8) * (1.0 / 16777216.0) - 0.5, + ((uint)hash.y >> 8) * (1.0 / 16777216.0) - 0.5, + ((uint)hash.z >> 8) * (1.0 / 16777216.0) - 0.5); + float3 v = local - gridStep - jitter; + float falloff = exp(-2.0 * dot(v, v)); + float w = max(0.0, falloff - 0.01111); + float3 dw = w > 0.0 ? (-4.0 * falloff) * v : float3(0.0, 0.0, 0.0); + float phase = dot(v, side) + phaseOffset; + float c = cos(phase); + float s = sin(phase); + + sumW += w; + sumCos += w * c; + sumSin += w * s; + gradW += dw; + gradCos += dw * c - (w * s) * side; + gradSin += dw * s + (w * c) * side; + } + } + } + + float invW = 1.0 / sumW; + float rawCos = sumCos * invW; + float rawSin = sumSin * invW; + float3 gradRawCos = (gradCos - rawCos * gradW) * invW; + float3 gradRawSin = (gradSin - rawSin * gradW) * invW; + float rawLength = sqrt(rawCos * rawCos + rawSin * rawSin); + float floorLength = 1.0 - normalization; + + CultPhasor result; + if (rawLength > floorLength) + { + float outCos = rawCos / rawLength; + float outSin = rawSin / rawLength; + float3 gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength; + result.cos = float4((gradRawCos - outCos * gradLength) / rawLength, outCos); + result.sin = float4((gradRawSin - outSin * gradLength) / rawLength, outSin); + } + else + { + result.cos = float4(gradRawCos / floorLength, rawCos / floorLength); + result.sin = float4(gradRawSin / floorLength, rawSin / floorLength); + } + + return result; +} + float cultmath_value_noise(float2 position) { float2 cell = floor(position); diff --git a/packages/cultmath/src/CultMath/CultPhasor.cs b/packages/cultmath/src/CultMath/CultPhasor.cs new file mode 100644 index 00000000..0a28faa8 --- /dev/null +++ b/packages/cultmath/src/CultMath/CultPhasor.cs @@ -0,0 +1,21 @@ +namespace CultMath; + +/// +/// The value-and-gradient result of +/// (invariant 8's struct return shape). Both fields use the float4(gradient.xyz, value.w) +/// layout of a bare value-and-gradient return. +/// +public struct CultPhasor +{ + /// (∇cos, cos): the normalized cosine wave and its gradient. + public float4 cos; + + /// (∇sin, sin): the normalized sine wave and its gradient. + public float4 sin; + + public CultPhasor(float4 cos, float4 sin) + { + this.cos = cos; + this.sin = sin; + } +} diff --git a/packages/cultmath/src/CultMath/CultPhasor.cs.meta b/packages/cultmath/src/CultMath/CultPhasor.cs.meta new file mode 100644 index 00000000..d586aa3f --- /dev/null +++ b/packages/cultmath/src/CultMath/CultPhasor.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 807ae2c8a758cf31896da12f1e5047bd diff --git a/packages/cultmath/src/CultMath/math.cs b/packages/cultmath/src/CultMath/math.cs index 23916b6b..430325ea 100644 --- a/packages/cultmath/src/CultMath/math.cs +++ b/packages/cultmath/src/CultMath/math.cs @@ -968,6 +968,104 @@ public static partial class math // dxc and C#, bit-exact on both sides. internal static float cellular_unit(int bits) => ((uint)bits >> 8) * (1.0f / 16777216.0f); + // Johansen's Phacelle noise (phase + cell), generalised to 3D cells with a caller-supplied + // stripe wave vector. Each of the 4x4x4 cells around p carries one jittered feature point + // (jitter in [-0.5, 0.5) per axis, from pcg3d over the integer cell coordinate, the same + // integer-hashing rule as cellular) and one cosine/sine wave that runs along `side` (the wave + // vector: its length is 2*pi times the stripe frequency, and the caller chooses which + // perpendicular it is; on a unit sphere that is cross(p_hat, flow)). The waves are blended with + // the bell weight max(0, exp(-2 d^2) - 0.01111) (exactly 0 at distance 1.5, the nearest a + // feature outside the 4x4x4 block can be in 2D, which is what keeps cell-grid lines out), and + // the blended phasor is normalized: its length is stretched to 1 wherever the raw length is at + // least 1 - normalization, and scaled by 1 / (1 - normalization) below that. `offset` is the + // phase in cycles (1 is a full turn). Value and normalization follow Johansen's + // PhacelleNoise as ported in lpmitchell/AdvancedTerrainErosion. + // + // The gradient is exact: it carries the weight derivative (grad w = -4 v exp(-2 |v|^2) where + // w > 0) and the normalization, not upstream's convention grad(cos) ~ -sin * side, which + // ignores the weight gradient (Asura probe: 16.8 deg off at the median, 97.8 at the 95th + // percentile; the exact form is 0.00 and 0.03). With P = sum(w * (cos, sin)) and W = sum(w), + // I = P / W has grad I = (grad P - I * grad W) / W; above the threshold the output is I / |I| + // with gradient (grad I - out * grad|I|) / |I|, below it grad I / (1 - normalization). The + // gradient is discontinuous where |I| crosses the threshold, and where a cell's weight leaves + // its support (distance 1.4998): there grad w jumps from -4 * 0.01111 * v (length about 0.067) to 0. The value is + // continuous at both. + // + // Pruning is exact: per axis, |local - g| - 0.5 (floored at 0) lower-bounds |v| for any jitter, + // so a cell whose summed square is at least 2.25 has weight exactly 0 (exp(-4.5) = 0.011109 < + // 0.01111) and contributes nothing to any sum. Only the internal overload can switch the prune + // off; PhacelleTests uses it to pin bit equality. About 14 of 64 cells survive on average. + public static CultPhasor phacelle(float3 p, float3 side, float offset, float normalization) => + phacelle(p, side, offset, normalization, true); + + internal static CultPhasor phacelle(float3 p, float3 side, float offset, float normalization, bool prune) + { + var cell = floor(p); + var local = p - cell; + var phaseOffset = offset * TAU; + var sumW = 0.0f; + var sumCos = 0.0f; + var sumSin = 0.0f; + var gradW = new float3(0.0f, 0.0f, 0.0f); + var gradCos = new float3(0.0f, 0.0f, 0.0f); + var gradSin = new float3(0.0f, 0.0f, 0.0f); + + for (var dz = -1; dz <= 2; dz++) + { + var gz = max(0.0f, abs(local.z - (float)dz) - 0.5f); + for (var dy = -1; dy <= 2; dy++) + { + var gy = max(0.0f, abs(local.y - (float)dy) - 0.5f); + for (var dx = -1; dx <= 2; dx++) + { + var gx = max(0.0f, abs(local.x - (float)dx) - 0.5f); + if (prune && gx * gx + gy * gy + gz * gz >= 2.25f) + continue; + + var gridStep = new float3(dx, dy, dz); + var hash = pcg3d(int3(cell + gridStep)); + var jitter = new float3(cellular_unit(hash.x) - 0.5f, cellular_unit(hash.y) - 0.5f, cellular_unit(hash.z) - 0.5f); + var v = local - gridStep - jitter; + var falloff = exp(-2.0f * dot(v, v)); + var w = max(0.0f, falloff - 0.01111f); + var dw = w > 0.0f ? (-4.0f * falloff) * v : new float3(0.0f, 0.0f, 0.0f); + var phase = dot(v, side) + phaseOffset; + var c = cos(phase); + var s = sin(phase); + + sumW += w; + sumCos += w * c; + sumSin += w * s; + gradW += dw; + gradCos += dw * c - (w * s) * side; + gradSin += dw * s + (w * c) * side; + } + } + } + + var invW = 1.0f / sumW; + var rawCos = sumCos * invW; + var rawSin = sumSin * invW; + var gradRawCos = (gradCos - rawCos * gradW) * invW; + var gradRawSin = (gradSin - rawSin * gradW) * invW; + var rawLength = sqrt(rawCos * rawCos + rawSin * rawSin); + var floorLength = 1.0f - normalization; + + if (rawLength > floorLength) + { + var outCos = rawCos / rawLength; + var outSin = rawSin / rawLength; + var gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength; + return new CultPhasor( + new float4((gradRawCos - outCos * gradLength) / rawLength, outCos), + new float4((gradRawSin - outSin * gradLength) / rawLength, outSin)); + } + + return new CultPhasor( + new float4(gradRawCos / floorLength, rawCos / floorLength), + new float4(gradRawSin / floorLength, rawSin / floorLength)); + } + public static float value_noise(float2 position) { var cell = floor(position); diff --git a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs index be8a3df4..ad690759 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs @@ -30,6 +30,7 @@ public sealed class HlslSourceCompatibilityTests Assert.NotNull(ShaderFunction(shader, "cultmath_snoise_grad", typeof(float3))); Assert.NotNull(ShaderFunction(shader, "cultmath_fbm_grad", typeof(float3), typeof(int), typeof(float), typeof(float))); Assert.NotNull(ShaderFunction(shader, "cultmath_ridged_grad", typeof(float3), typeof(int), typeof(float), typeof(float))); + Assert.NotNull(ShaderFunction(shader, "cultmath_phacelle", typeof(float3), typeof(float3), typeof(float), typeof(float))); } // HLSL functions carry no access modifier, so they compile as private instance methods. @@ -130,7 +131,7 @@ public sealed class HlslSourceCompatibilityTests } Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}"); - Assert.Equal(33, compared.Count); + Assert.Equal(34, compared.Count); } /// diff --git a/packages/cultmath/tests/CultMath.Tests/PhacelleTests.cs b/packages/cultmath/tests/CultMath.Tests/PhacelleTests.cs new file mode 100644 index 00000000..15d1a89a --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/PhacelleTests.cs @@ -0,0 +1,275 @@ +using System; +using CultMath; +using Xunit; + +namespace CultMath.Tests; + +public sealed class PhacelleTests +{ + // Every tolerance below is at least 5x the largest error the same test measures at its seed. + private const float GradientEps = 1.0e-3f; + private const float GradientTolerance = 5.0e-3f; + private const double ValueTolerance = 1.6e-5; + private const double ThresholdBand = 0.02; + + private static float3 RandomPoint(System.Random random, float extent = 3.0f) => new( + (random.NextSingle() * 2.0f - 1.0f) * extent, (random.NextSingle() * 2.0f - 1.0f) * extent, (random.NextSingle() * 2.0f - 1.0f) * extent); + + // A wave vector of length 2*pi times a stripe frequency in [0.5, 1.5], in a random direction. + private static float3 RandomSide(System.Random random) + { + float3 direction; + do + { + direction = new float3(random.NextSingle() - 0.5f, random.NextSingle() - 0.5f, random.NextSingle() - 0.5f); + } + while (math.length(direction) < 0.1f); + + return math.normalize(direction) * (math.TAU * (0.5f + random.NextSingle())); + } + + // ---- An independent double-precision evaluation: all 64 cells, no prune, no gradient. ---- + // + // It shares only the hash (pcg3d over the integer cell coordinate, top 24 bits) and the constants + // the function is specified by: weight max(0, exp(-2 d^2) - 0.01111), jitter in [-0.5, 0.5). + + private readonly record struct Reference(double Cos, double Sin, double RawLength, double SupportEdgeGap); + + private static Reference Evaluate(float3 p, float3 side, float offset, float normalization) + { + double fx = Math.Floor(p.x), fy = Math.Floor(p.y), fz = Math.Floor(p.z); + double lx = p.x - fx, ly = p.y - fy, lz = p.z - fz; + double sumW = 0.0, sumCos = 0.0, sumSin = 0.0, edgeGap = double.MaxValue; + var supportRadius = Math.Sqrt(Math.Log(1.0 / 0.01111) / 2.0); + for (var dz = -1; dz <= 2; dz++) + for (var dy = -1; dy <= 2; dy++) + for (var dx = -1; dx <= 2; dx++) + { + var hash = math.pcg3d(new int3((int)fx + dx, (int)fy + dy, (int)fz + dz)); + double vx = lx - dx - (math.cellular_unit(hash.x) - 0.5); + double vy = ly - dy - (math.cellular_unit(hash.y) - 0.5); + double vz = lz - dz - (math.cellular_unit(hash.z) - 0.5); + var d2 = vx * vx + vy * vy + vz * vz; + edgeGap = Math.Min(edgeGap, Math.Abs(Math.Sqrt(d2) - supportRadius)); + var w = Math.Max(0.0, Math.Exp(-2.0 * d2) - 0.01111); + var phase = vx * side.x + vy * side.y + vz * side.z + (double)offset * math.TAU; + sumW += w; + sumCos += w * Math.Cos(phase); + sumSin += w * Math.Sin(phase); + } + + double rawCos = sumCos / sumW, rawSin = sumSin / sumW; + var raw = Math.Sqrt(rawCos * rawCos + rawSin * rawSin); + var length = Math.Max(1.0 - normalization, raw); + return new Reference(rawCos / length, rawSin / length, raw, edgeGap); + } + + private static float3 CentralDifference(Func f, float3 p, float eps) => new( + (f(p + new float3(eps, 0.0f, 0.0f)) - f(p - new float3(eps, 0.0f, 0.0f))) / (2.0f * eps), + (f(p + new float3(0.0f, eps, 0.0f)) - f(p - new float3(0.0f, eps, 0.0f))) / (2.0f * eps), + (f(p + new float3(0.0f, 0.0f, eps)) - f(p - new float3(0.0f, 0.0f, eps))) / (2.0f * eps)); + + private static float MaxComponentError(float3 numeric, float4 analytic) => MathF.Max( + MathF.Abs(numeric.x - analytic.x), MathF.Max(MathF.Abs(numeric.y - analytic.y), MathF.Abs(numeric.z - analytic.z))); + + // ---- Value: the function is Johansen's, on 3D cells ---- + + [Theory] + [InlineData(0.0f)] + [InlineData(0.5f)] + [InlineData(0.9f)] + [InlineData(1.0f)] + public void ValueMatchesTheIndependentReferenceAndNormalizesAboveTheThreshold(float normalization) + { + var random = new System.Random(0x7ACE1 + (int)(normalization * 10.0f)); + var maxError = 0.0; + var above = 0; + var below = 0; + for (var i = 0; i < 3000; i++) + { + var p = RandomPoint(random, 40.0f); + var side = RandomSide(random); + var offset = random.NextSingle(); + var actual = math.phacelle(p, side, offset, normalization); + var expected = Evaluate(p, side, offset, normalization); + + maxError = Math.Max(maxError, Math.Max(Math.Abs(actual.cos.w - expected.Cos), Math.Abs(actual.sin.w - expected.Sin))); + if (expected.RawLength >= 1.0 - normalization + ThresholdBand) + { + // The output magnitude is 1 wherever the raw magnitude reaches 1 - normalization. + above++; + var magnitude = Math.Sqrt((double)actual.cos.w * actual.cos.w + (double)actual.sin.w * actual.sin.w); + Assert.True(Math.Abs(magnitude - 1.0) < ValueTolerance, $"magnitude {magnitude} above the threshold at {p}"); + } + else if (expected.RawLength <= 1.0 - normalization - ThresholdBand) + { + below++; + } + } + + Assert.True(maxError < ValueTolerance, $"max value error against the reference {maxError}"); + if (normalization > 0.0f) + Assert.True(above >= 500, $"only {above} points above the threshold"); + if (normalization is > 0.0f and < 1.0f) + Assert.True(below >= 50, $"only {below} points below the threshold"); + } + + [Fact] + public void ValueIsPeriodicInOffsetAndAQuarterCycleRotatesCosIntoSin() + { + // offset is in cycles: a quarter turn maps cos to -sin and sin to cos, and a whole turn is the + // identity. The weights and the normalization do not depend on the phase, so this is exact + // up to float rounding. + var random = new System.Random(0x7ACE2); + for (var i = 0; i < 200; i++) + { + var p = RandomPoint(random); + var side = RandomSide(random); + var offset = random.NextSingle(); + var basePhasor = math.phacelle(p, side, offset, 0.5f); + var quarter = math.phacelle(p, side, offset + 0.25f, 0.5f); + var whole = math.phacelle(p, side, offset + 1.0f, 0.5f); + + Assert.True(MathF.Abs(quarter.cos.w + basePhasor.sin.w) < 2.0e-5f); + Assert.True(MathF.Abs(quarter.sin.w - basePhasor.cos.w) < 2.0e-5f); + Assert.True(MathF.Abs(whole.cos.w - basePhasor.cos.w) < 2.0e-5f); + Assert.True(MathF.Abs(whole.sin.w - basePhasor.sin.w) < 2.0e-5f); + } + } + + // ---- Invariant 8: the gradient is exact ---- + + [Fact] + public void GradientOfCosAndSinMatchesCentralDifferencesOnBothSidesOfTheNormalizationThreshold() + { + // Above the threshold the output is I / |I|; below it, I / (1 - normalization). The gradient + // is discontinuous where |I| crosses the threshold, and at the edge of each cell's weight support + // (distance 1.4998, where max(0, ...) starts: a jump of about 4 d 0.01111 in grad w), so a band + // around each is excluded (by the + // reference's raw length, since the output alone cannot say how close it is). + const float normalization = 0.5f; + var random = new System.Random(0x7ACE3); + var above = 0; + var below = 0; + var maxError = 0.0f; + for (var attempts = 0; attempts < 20000 && (above < 150 || below < 150); attempts++) + { + var p = RandomPoint(random); + var side = RandomSide(random); + var offset = random.NextSingle(); + var reference = Evaluate(p, side, offset, normalization); + var distance = reference.RawLength - (1.0 - normalization); + if (Math.Abs(distance) < ThresholdBand || reference.SupportEdgeGap < 3.0 * GradientEps) + continue; + if (distance > 0.0 ? above >= 150 : below >= 150) + continue; + if (distance > 0.0) above++; else below++; + + var analytic = math.phacelle(p, side, offset, normalization); + var numericCos = CentralDifference(q => math.phacelle(q, side, offset, normalization).cos.w, p, GradientEps); + var numericSin = CentralDifference(q => math.phacelle(q, side, offset, normalization).sin.w, p, GradientEps); + maxError = MathF.Max(maxError, MathF.Max(MaxComponentError(numericCos, analytic.cos), MaxComponentError(numericSin, analytic.sin))); + } + + Assert.True(above >= 150 && below >= 150, $"{above} points above and {below} below the threshold"); + Assert.True(maxError < GradientTolerance, $"max gradient error against central differences {maxError}"); + } + + // ---- Stripes run along the flow ---- + + [Fact] + public void OnTheUnitSphereStripesRunAlongTheFlow() + { + // side = cross(p_hat, flow) * 2*pi is tangent to the sphere and perpendicular to the flow, so the + // cosine should change across the stripes (along side) and hardly at all along them (along + // flow). Asura's probe measured 0.212 for this ratio, against 0.202 for the planar upstream. + const float radiusInCells = 10.0f; + var random = new System.Random(0x7ACE4); + var reference = math.normalize(new float3(0.3f, 0.5f, 0.8f)); + var alongFlow = 0.0; + var acrossFlow = 0.0; + var used = 0; + while (used < 2000) + { + var hat = math.normalize(new float3(random.NextSingle() - 0.5f, random.NextSingle() - 0.5f, random.NextSingle() - 0.5f)); + var tangent = reference - hat * math.dot(hat, reference); + if (math.length(tangent) < 0.2f) + continue; + + var flow = math.normalize(tangent); + var sideDirection = math.cross(hat, flow); + var phasor = math.phacelle(hat * radiusInCells, sideDirection * math.TAU, 0.25f, 0.5f); + var gradient = new float3(phasor.cos.x, phasor.cos.y, phasor.cos.z); + alongFlow += Math.Abs(math.dot(gradient, flow)); + acrossFlow += Math.Abs(math.dot(gradient, sideDirection)); + used++; + } + + var ratio = alongFlow / acrossFlow; + Assert.True(ratio <= 0.25, $"mean|d cos/d flow| / mean|d cos/d side| = {ratio}"); + } + + // ---- Continuity ---- + + [Fact] + public void ValueStepDifferenceShrinksTenfoldPerTenfoldStepIncludingAcrossCellPlanes() + { + // Pairs straddle an integer plane on each axis (where the 4x4x4 window moves by one cell and a + // weight that did not vanish at the window edge would jump) and also start at random points. + // Lipschitz continuity means the largest difference over a step scales with the step. + var random = new System.Random(0x7ACE5); + var steps = new[] { 1.0e-3f, 1.0e-4f, 1.0e-5f }; + var maxDifference = new float[steps.Length]; + for (var i = 0; i < 3000; i++) + { + var p = RandomPoint(random); + var side = RandomSide(random); + var offset = random.NextSingle(); + var axis = i % 3; + if (i % 2 == 0) + p = SetAxis(p, axis, MathF.Round(GetAxis(p, axis))); + for (var k = 0; k < steps.Length; k++) + { + var low = math.phacelle(SetAxis(p, axis, GetAxis(p, axis) - steps[k] * 0.5f), side, offset, 0.5f); + var high = math.phacelle(SetAxis(p, axis, GetAxis(p, axis) + steps[k] * 0.5f), side, offset, 0.5f); + maxDifference[k] = MathF.Max(maxDifference[k], MathF.Max(MathF.Abs(low.cos.w - high.cos.w), MathF.Abs(low.sin.w - high.sin.w))); + } + } + + Assert.True(maxDifference[1] * 8.0f <= maxDifference[0], $"1e-3: {maxDifference[0]}, 1e-4: {maxDifference[1]}"); + Assert.True(maxDifference[2] * 8.0f <= maxDifference[1], $"1e-4: {maxDifference[1]}, 1e-5: {maxDifference[2]}"); + } + + private static float GetAxis(float3 v, int axis) => axis == 0 ? v.x : axis == 1 ? v.y : v.z; + + private static float3 SetAxis(float3 v, int axis, float value) => + axis == 0 ? new float3(value, v.y, v.z) : axis == 1 ? new float3(v.x, value, v.z) : new float3(v.x, v.y, value); + + // ---- Exact pruning ---- + + [Fact] + public void PrunedOutputIsBitEqualToTheUnprunedSum() + { + var random = new System.Random(0x7ACE6); + for (var i = 0; i < 20000; i++) + { + var p = RandomPoint(random, 50.0f); + var side = RandomSide(random); + var offset = random.NextSingle(); + var normalization = random.NextSingle(); + var pruned = math.phacelle(p, side, offset, normalization); + var unpruned = math.phacelle(p, side, offset, normalization, false); + + AssertBitEqual(unpruned.cos, pruned.cos); + AssertBitEqual(unpruned.sin, pruned.sin); + } + } + + private static void AssertBitEqual(float4 expected, float4 actual) + { + Assert.Equal(BitConverter.SingleToInt32Bits(expected.x), BitConverter.SingleToInt32Bits(actual.x)); + Assert.Equal(BitConverter.SingleToInt32Bits(expected.y), BitConverter.SingleToInt32Bits(actual.y)); + Assert.Equal(BitConverter.SingleToInt32Bits(expected.z), BitConverter.SingleToInt32Bits(actual.z)); + Assert.Equal(BitConverter.SingleToInt32Bits(expected.w), BitConverter.SingleToInt32Bits(actual.w)); + } +} diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl index 50c5a746..d0dd22ae 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl @@ -503,6 +503,92 @@ CultCellular cultmath_cellular(float3 p) return result; } +// Invariant 8's struct return shape for cultmath_phacelle, compared field by field and bit for bit +// against CultMath.CultPhasor, same field order. +struct CultPhasor +{ + float4 cos; + float4 sin; +}; + +// Johansen's Phacelle noise on 3D cells. See math.phacelle's comment (math.cs) for the stripe wave +// vector, the weight, the normalization, the exact gradient and why the prune is exact. Mirrors +// the C# body operation for operation; the C#-only internal overload that switches the prune off +// has no shader twin. +CultPhasor cultmath_phacelle(float3 p, float3 side, float offset, float normalization) +{ + float3 cell = floor(p); + float3 local = p - cell; + float phaseOffset = offset * CULTMATH_TAU; + float sumW = 0.0; + float sumCos = 0.0; + float sumSin = 0.0; + float3 gradW = float3(0.0, 0.0, 0.0); + float3 gradCos = float3(0.0, 0.0, 0.0); + float3 gradSin = float3(0.0, 0.0, 0.0); + + for (int dz = -1; dz <= 2; dz++) + { + float gz = max(0.0, abs(local.z - (float)dz) - 0.5); + for (int dy = -1; dy <= 2; dy++) + { + float gy = max(0.0, abs(local.y - (float)dy) - 0.5); + for (int dx = -1; dx <= 2; dx++) + { + float gx = max(0.0, abs(local.x - (float)dx) - 0.5); + if (gx * gx + gy * gy + gz * gz >= 2.25) + continue; + + float3 gridStep = float3(dx, dy, dz); + int3 hash = cultmath_pcg3d(int3(cell + gridStep)); + float3 jitter = float3( + ((uint)hash.x >> 8) * (1.0 / 16777216.0) - 0.5, + ((uint)hash.y >> 8) * (1.0 / 16777216.0) - 0.5, + ((uint)hash.z >> 8) * (1.0 / 16777216.0) - 0.5); + float3 v = local - gridStep - jitter; + float falloff = exp(-2.0 * dot(v, v)); + float w = max(0.0, falloff - 0.01111); + float3 dw = w > 0.0 ? (-4.0 * falloff) * v : float3(0.0, 0.0, 0.0); + float phase = dot(v, side) + phaseOffset; + float c = cos(phase); + float s = sin(phase); + + sumW += w; + sumCos += w * c; + sumSin += w * s; + gradW += dw; + gradCos += dw * c - (w * s) * side; + gradSin += dw * s + (w * c) * side; + } + } + } + + float invW = 1.0 / sumW; + float rawCos = sumCos * invW; + float rawSin = sumSin * invW; + float3 gradRawCos = (gradCos - rawCos * gradW) * invW; + float3 gradRawSin = (gradSin - rawSin * gradW) * invW; + float rawLength = sqrt(rawCos * rawCos + rawSin * rawSin); + float floorLength = 1.0 - normalization; + + CultPhasor result; + if (rawLength > floorLength) + { + float outCos = rawCos / rawLength; + float outSin = rawSin / rawLength; + float3 gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength; + result.cos = float4((gradRawCos - outCos * gradLength) / rawLength, outCos); + result.sin = float4((gradRawSin - outSin * gradLength) / rawLength, outSin); + } + else + { + result.cos = float4(gradRawCos / floorLength, rawCos / floorLength); + result.sin = float4(gradRawSin / floorLength, rawSin / floorLength); + } + + return result; +} + float cultmath_value_noise(float2 position) { float2 cell = floor(position); From 9f18a432cba499a42da96f3fa1b3a465c4b7fdc9 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 23:51:04 +0200 Subject: [PATCH 2/4] CultMath: document phacelle and record Johansen's MPL-2.0 provenance design.md names CultPhasor as the second struct-return consumer and describes phacelle (stripe wave vector, weight, normalization, exact gradient and its two discontinuity classes, exact pruning). THIRD-PARTY-NOTICES.md, in both the package root and the Unity package, cites Rune Skovbo Johansen's Phacelle Noise (MPL-2.0) and records that the 3D-cell generalisation, the caller-supplied side vector, the exact gradient and the exact pruning are new. Co-Authored-By: Claude Opus 5.5 --- packages/cultmath/THIRD-PARTY-NOTICES.md | 15 ++++++++++++ packages/cultmath/docs/design.md | 23 +++++++++++++++++-- .../THIRD-PARTY-NOTICES.md | 15 ++++++++++++ 3 files changed, 51 insertions(+), 2 deletions(-) diff --git a/packages/cultmath/THIRD-PARTY-NOTICES.md b/packages/cultmath/THIRD-PARTY-NOTICES.md index 1de744a7..73f2461a 100644 --- a/packages/cultmath/THIRD-PARTY-NOTICES.md +++ b/packages/cultmath/THIRD-PARTY-NOTICES.md @@ -19,3 +19,18 @@ permutation and scale as `snoise(float3)`. - Copyright 2011 Ashima Arts. - Upstream: - License: MIT. + +## Phacelle Noise + +`math.phacelle` and `cultmath_phacelle` follow Rune Skovbo Johansen's Phacelle +noise (Shadertoy `t3dyWl`; erosion filter `wXcfWn`; blog, March 2026), read +from the C# port `PhacelleNoise` in `lpmitchell/AdvancedTerrainErosion` +(`AdvancedTerrainErosion.cs`, converted to burstable C# by Luke Mitchell, +2026). The cell weight, the phase blend and the normalization follow it. The +3D-cell generalisation, the caller-supplied stripe wave vector, the exact +gradient and the exact pruning are new here. + +- Copyright (c) 2025 Rune Skovbo Johansen. +- Upstream: +- License: MPL-2.0 (). CultMath's Unity package + is MPL-2.0, so the licence is compatible. diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index e14927fe..0858ee8c 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -128,8 +128,8 @@ whose HLSL return type is a plain struct of `float`/`floatN`/`int`/`intN` fields is compared to its C# counterpart field by field, recursing into each field's own components, down to bit-for-bit scalars, matching each leaf's declaration path (field names, not just position) and declared type as well -as its value. `CultCellular` (`math.cs`, `cellular`) is the first and, for -now, only consumer. A struct return with a mismatched, reordered-by-name, or +as its value. `CultCellular` (`math.cs`, `cellular`) and `CultPhasor` +(`math.cs`, `phacelle`) are its consumers. A struct return with a mismatched, reordered-by-name, or retyped field does not get a second comparison path or a shape of its own; it fails the same walk that already handles a bare vector return. @@ -264,6 +264,25 @@ together, laid out as `float4(gradient.xyz, value.w)`, with no value-only twin 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 these two mirrors. +- `phacelle(float3 p, float3 side, float offset, float normalization)` is Rune + Skovbo Johansen's Phacelle noise ("phase + cell": a stripe pattern built by + blending one cosine/sine wave per jittered cell), generalised to 3D cells. It + returns `CultPhasor { cos, sin }`, each `(∇, value)`. `side` is the stripe + wave vector, whose length is 2π times the stripe frequency; in 3D the + perpendicular to the flow is not unique, so the caller chooses it (on a unit + sphere, `cross(p̂, flow)`). `offset` is the phase in cycles. It visits 4×4×4 + cells with `pcg3d`-hashed jitter in `[-0.5, 0.5)`, weights each wave by + `max(0, exp(-2d²) - 0.01111)`, and normalizes the blended phasor as upstream + does (length 1 where the raw length is at least `1 - normalization`, scaled + by `1 / (1 - normalization)` below). The gradient is exact, weight + derivatives and normalization included, rather than upstream's + `∇cos ≈ -sin·side`, which ignores the weight gradient. It is discontinuous + where the raw length crosses the normalization threshold and where a cell's + weight leaves its support (distance 1.4998); the value is continuous at both. + Cells whose per-axis lower bound `Σ max(0, |local - g| - 0.5)²` is at least + 2.25 have weight exactly 0 and are skipped, which changes no output bit; + `PhacelleTests` pins that against an internal overload with the prune + disabled. No 2D variant exists; nothing consumes one. Integer hashing uses the PCG hashes from Jarzynski and Olano, "Hash Functions for GPU Rendering" (JCGT 9(3), 2020): `pcg(uint)` is O'Neill's RXS-M-XS 32/32 diff --git a/packages/cultmath/unity/org.gamecult.cultmath/THIRD-PARTY-NOTICES.md b/packages/cultmath/unity/org.gamecult.cultmath/THIRD-PARTY-NOTICES.md index 1de744a7..73f2461a 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/THIRD-PARTY-NOTICES.md +++ b/packages/cultmath/unity/org.gamecult.cultmath/THIRD-PARTY-NOTICES.md @@ -19,3 +19,18 @@ permutation and scale as `snoise(float3)`. - Copyright 2011 Ashima Arts. - Upstream: - License: MIT. + +## Phacelle Noise + +`math.phacelle` and `cultmath_phacelle` follow Rune Skovbo Johansen's Phacelle +noise (Shadertoy `t3dyWl`; erosion filter `wXcfWn`; blog, March 2026), read +from the C# port `PhacelleNoise` in `lpmitchell/AdvancedTerrainErosion` +(`AdvancedTerrainErosion.cs`, converted to burstable C# by Luke Mitchell, +2026). The cell weight, the phase blend and the normalization follow it. The +3D-cell generalisation, the caller-supplied stripe wave vector, the exact +gradient and the exact pruning are new here. + +- Copyright (c) 2025 Rune Skovbo Johansen. +- Upstream: +- License: MPL-2.0 (). CultMath's Unity package + is MPL-2.0, so the licence is compatible. From 2a3cab1bd7ab5f2532f8a2ff32f328d057ae2de5 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 23:58:43 +0200 Subject: [PATCH 3/4] CultMath tests: a dense sweep pins cultmath_phacelle against C# bit for bit The generic mirror comparison passes about a hundred inputs through the shader and missed a mirror whose prune cut was 2.0 instead of 2.25, because only a rare cell has a nonzero weight under that cut. Twenty thousand random arguments compare the shader text and math.phacelle field by field and kill that mutant. Co-Authored-By: Claude Opus 5.5 --- .../HlslSourceCompatibilityTests.cs | 29 +++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs index ad690759..2cf737df 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs @@ -134,6 +134,35 @@ public sealed class HlslSourceCompatibilityTests Assert.Equal(34, compared.Count); } + /// + /// The generic comparison above passes about a hundred inputs through cultmath_phacelle, and + /// its prune only changes an output for the rare cell whose weight is nonzero but whose lower + /// bound sits just under the 2.25 cut (Soul's F6 lesson for cellular): a mirror whose cut is 2.0 + /// survives a hundred inputs. A dense sweep of realistic arguments closes that gap. + /// + [Fact] + public void PhacelleMirrorMatchesCSharpBitForBitAcrossADenseSweep() + { + var (assembly, errors) = CompileShaderMirror(); + Assert.True(assembly is not null, string.Join(Environment.NewLine, errors)); + var shaderType = assembly!.GetType("CultMathHlsl.HlslShader")!; + var shader = Activator.CreateInstance(shaderType); + var mirror = ShaderFunction(shaderType, "cultmath_phacelle", typeof(float3), typeof(float3), typeof(float), typeof(float))!; + + var random = new System.Random(0x5EED2); + float Next(float extent) => (random.NextSingle() * 2.0f - 1.0f) * extent; + for (var i = 0; i < 20000; i++) + { + var p = new float3(Next(50.0f), Next(50.0f), Next(50.0f)); + var side = new float3(Next(9.0f), Next(9.0f), Next(9.0f)); + var offset = random.NextSingle(); + var normalization = random.NextSingle(); + var expected = Values(phacelle(p, side, offset, normalization)); + var actual = Values(mirror.Invoke(shader, new object[] { p, side, offset, normalization })!); + Assert.True(BitwiseEqual(expected, actual), $"mirror differs from C# at p = {p}, side = {side}"); + } + } + /// /// A hand-built pair of structs standing in for a mirror struct return and its C# counterpart /// (invariant 8's one struct return shape, e.g. ): same field shape, From 0acc3728ed521a847d589f8f0cd9e9202f64ca96 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Wed, 30 Sep 2026 00:34:27 +0200 Subject: [PATCH 4/4] CultMath: Phacelle port moves into MPL-2.0 files; prune cut pinned by boundary points MPL-2.0 is file-level and the Phacelle port sat inside MIT files. phacelle and CultPhasor now live in math.Phacelle.cs and CultPhasor.cs, and cultmath_phacelle in CultMath.Phacelle.hlsl (included by CultMath.hlsl, both tracked copies), each with the MPL header and Johansen / Luke Mitchell attribution. math stays one partial class, so math.phacelle is called exactly as before. The GameCult.Math package declares MIT AND MPL-2.0 and packs THIRD-PARTY-NOTICES.md and the new include; the Unity LICENSE.md points at CultLib; build-unity-package.ps1 requires the new include. The random sweeps never met a cell whose box bound sits just under the 2.25 cut while its weight is still nonzero, so cuts of 2.245, 2.24 and 2.2 survived. Six hard-coded points (found by an offline search) sit next to such cells; the C# pruned-vs-unpruned test and the HLSL mirror test both use them, with a witness check that re-derives the property from the hash. Docs: the prune keeps about 45 of 64 cells (14 is the nonzero-weight count), the support radius is 1.49998 (d^2 = 2.24995), normalization's range and saturation, W = 0 only for adversarial hashes, and the support-edge gradient seam size. Co-Authored-By: Claude Sonnet 5.5 --- packages/cultmath/README.md | 2 + packages/cultmath/THIRD-PARTY-NOTICES.md | 9 +- packages/cultmath/docs/design.md | 26 +++- .../cultmath/scripts/build-unity-package.ps1 | 4 +- .../cultmath/shaders/CultMath.Phacelle.hlsl | 99 ++++++++++++++ .../shaders/CultMath.Phacelle.hlsl.meta | 7 + packages/cultmath/shaders/CultMath.hlsl | 86 +----------- .../cultmath/src/CultMath/CultMath.csproj | 4 + packages/cultmath/src/CultMath/CultPhasor.cs | 6 + .../cultmath/src/CultMath/math.Phacelle.cs | 127 ++++++++++++++++++ .../src/CultMath/math.Phacelle.cs.meta | 2 + packages/cultmath/src/CultMath/math.cs | 98 -------------- .../tests/CultMath.Tests/HlslMirrorTests.cs | 15 ++- .../HlslSourceCompatibilityTests.cs | 41 +++++- .../tests/CultMath.Tests/PhacelleTests.cs | 67 ++++++++- .../unity/org.gamecult.cultmath/LICENSE.md | 4 +- .../unity/org.gamecult.cultmath/README.md | 2 +- .../Shaders/CultMath.Phacelle.hlsl | 99 ++++++++++++++ .../Shaders/CultMath.Phacelle.hlsl.meta | 7 + .../Shaders/CultMath.hlsl | 86 +----------- .../THIRD-PARTY-NOTICES.md | 9 +- 21 files changed, 515 insertions(+), 285 deletions(-) create mode 100644 packages/cultmath/shaders/CultMath.Phacelle.hlsl create mode 100644 packages/cultmath/shaders/CultMath.Phacelle.hlsl.meta create mode 100644 packages/cultmath/src/CultMath/math.Phacelle.cs create mode 100644 packages/cultmath/src/CultMath/math.Phacelle.cs.meta create mode 100644 packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Phacelle.hlsl create mode 100644 packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Phacelle.hlsl.meta diff --git a/packages/cultmath/README.md b/packages/cultmath/README.md index 13f17973..f70df445 100644 --- a/packages/cultmath/README.md +++ b/packages/cultmath/README.md @@ -33,6 +33,8 @@ shader semantics. - `shaders/CultMath.hlsl`, a canonical HLSL mirror include for shader-side parity. HLSL already owns `float2`, `float3`, and `float4`; the include exposes `cultmath_*` functions for shared semantics HLSL does not provide. + It includes `shaders/CultMath.Phacelle.hlsl` (MPL-2.0) from the same + directory, so ship both files. - `Voronoi.SampleTones`, a C# batch surface that calls the Rust `cultmath-core` native kernel when `cultmath_core` is available and falls back to the managed parity path otherwise. diff --git a/packages/cultmath/THIRD-PARTY-NOTICES.md b/packages/cultmath/THIRD-PARTY-NOTICES.md index 73f2461a..625d63a9 100644 --- a/packages/cultmath/THIRD-PARTY-NOTICES.md +++ b/packages/cultmath/THIRD-PARTY-NOTICES.md @@ -32,5 +32,10 @@ gradient and the exact pruning are new here. - Copyright (c) 2025 Rune Skovbo Johansen. - Upstream: -- License: MPL-2.0 (). CultMath's Unity package - is MPL-2.0, so the licence is compatible. +- License: MPL-2.0 (). MPL-2.0 is file-level, so it + covers exactly the files that carry the port and nothing else: + `src/CultMath/math.Phacelle.cs`, `src/CultMath/CultPhasor.cs`, + `shaders/CultMath.Phacelle.hlsl` and the Unity copy of that include. Each + carries the MPL-2.0 header and this attribution. The rest of CultMath is MIT, + so the `GameCult.Math` NuGet package declares `MIT AND MPL-2.0` and ships this + file. The Unity package as a whole is MPL-2.0. diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index 0858ee8c..0de8eb14 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -129,7 +129,7 @@ fields is compared to its C# counterpart field by field, recursing into each field's own components, down to bit-for-bit scalars, matching each leaf's declaration path (field names, not just position) and declared type as well as its value. `CultCellular` (`math.cs`, `cellular`) and `CultPhasor` -(`math.cs`, `phacelle`) are its consumers. A struct return with a mismatched, reordered-by-name, or +(`math.Phacelle.cs`, `phacelle`) are its consumers. A struct return with a mismatched, reordered-by-name, or retyped field does not get a second comparison path or a shape of its own; it fails the same walk that already handles a bare vector return. @@ -278,11 +278,27 @@ together, laid out as `float4(gradient.xyz, value.w)`, with no value-only twin derivatives and normalization included, rather than upstream's `∇cos ≈ -sin·side`, which ignores the weight gradient. It is discontinuous where the raw length crosses the normalization threshold and where a cell's - weight leaves its support (distance 1.4998); the value is continuous at both. + weight leaves its support (d² = 2.24995, distance 1.49998); the value is + continuous at both. At the support edge `∇w` jumps from `-4·0.01111·v` to 0, + a per-cell jump of length `0.0444·|v|` (0.0667 at `|v| = 1.5`); measured over + random points, the output gradient jumps by a median 1.5% of its length and + at most about 27%. `normalization` is meant to lie in `[0, 1]`: 1 or more + saturates (every output is stretched to length 1), a negative value scales + every output down by `1 / (1 - normalization)`, and the normalized gradient + scales as `1 / |I|` for the blended phasor `I`, so it grows where the blend + nearly cancels. The total weight `W` is positive for every `pcg3d` input (the + smallest found by hill-climbing is 0.056); `W = 0`, and a NaN output, needs + adversarial jitter no hash yields. Cells whose per-axis lower bound `Σ max(0, |local - g| - 0.5)²` is at least - 2.25 have weight exactly 0 and are skipped, which changes no output bit; - `PhacelleTests` pins that against an internal overload with the prune - disabled. No 2D variant exists; nothing consumes one. + 2.25 lie outside the support and have weight exactly 0, so they are skipped, + which changes no output bit. Cells whose bound sits just under 2.25 can still + carry a nonzero weight, so the cut cannot be lowered; `PhacelleTests` pins + the cut with points next to such cells, against an internal overload with the + prune disabled, and the HLSL mirror test pins it on the shader side. The prune + keeps about 45 of the 64 cells on average (a saving of about 30%); about 14 + have a nonzero weight. The port sits in its own MPL-2.0 files + (`math.Phacelle.cs`, `CultPhasor.cs`, `CultMath.Phacelle.hlsl`); see + `THIRD-PARTY-NOTICES.md`. No 2D variant exists; nothing consumes one. Integer hashing uses the PCG hashes from Jarzynski and Olano, "Hash Functions for GPU Rendering" (JCGT 9(3), 2020): `pcg(uint)` is O'Neill's RXS-M-XS 32/32 diff --git a/packages/cultmath/scripts/build-unity-package.ps1 b/packages/cultmath/scripts/build-unity-package.ps1 index 39e55813..228d735c 100644 --- a/packages/cultmath/scripts/build-unity-package.ps1 +++ b/packages/cultmath/scripts/build-unity-package.ps1 @@ -78,7 +78,9 @@ $requiredTemplateFiles = @( "Runtime\UnityBridge\UnityConversions.cs", "Runtime\UnityBridge\UnityConversions.cs.meta", "Shaders\CultMath.hlsl", - "Shaders\CultMath.hlsl.meta" + "Shaders\CultMath.hlsl.meta", + "Shaders\CultMath.Phacelle.hlsl", + "Shaders\CultMath.Phacelle.hlsl.meta" ) foreach ($relativePath in $requiredTemplateFiles) { if (-not (Test-Path -LiteralPath (Join-Path $templateRoot $relativePath) -PathType Leaf)) { diff --git a/packages/cultmath/shaders/CultMath.Phacelle.hlsl b/packages/cultmath/shaders/CultMath.Phacelle.hlsl new file mode 100644 index 00000000..1db0a2fc --- /dev/null +++ b/packages/cultmath/shaders/CultMath.Phacelle.hlsl @@ -0,0 +1,99 @@ +#ifndef CULTMATH_PHACELLE_HLSL +#define CULTMATH_PHACELLE_HLSL + +// This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of +// the MPL was not distributed with this file, You can obtain one at https://mozilla.org/MPL/2.0/. +// +// Phacelle noise. Copyright (c) 2025 Rune Skovbo Johansen (Shadertoy t3dyWl; MPL-2.0), read from the +// C# port PhacelleNoise in lpmitchell/AdvancedTerrainErosion, converted to burstable C# by Luke +// Mitchell (2026). Modified by GameCult: 3D cells, a caller-supplied stripe wave vector, an exact +// gradient and exact pruning. This file is an MPL-2.0 file; CultMath.hlsl, which includes it after +// its own definitions of CULTMATH_TAU and cultmath_pcg3d, is MIT (see THIRD-PARTY-NOTICES.md). + +// Invariant 8's struct return shape for cultmath_phacelle, compared field by field and bit for bit +// against CultMath.CultPhasor, same field order. +struct CultPhasor +{ + float4 cos; + float4 sin; +}; + +// Johansen's Phacelle noise on 3D cells. See math.phacelle's comment (math.Phacelle.cs) for the stripe wave +// vector, the weight, the normalization, the exact gradient and why the prune is exact. Mirrors +// the C# body operation for operation; the C#-only internal overload that switches the prune off +// has no shader twin. +CultPhasor cultmath_phacelle(float3 p, float3 side, float offset, float normalization) +{ + float3 cell = floor(p); + float3 local = p - cell; + float phaseOffset = offset * CULTMATH_TAU; + float sumW = 0.0; + float sumCos = 0.0; + float sumSin = 0.0; + float3 gradW = float3(0.0, 0.0, 0.0); + float3 gradCos = float3(0.0, 0.0, 0.0); + float3 gradSin = float3(0.0, 0.0, 0.0); + + for (int dz = -1; dz <= 2; dz++) + { + float gz = max(0.0, abs(local.z - (float)dz) - 0.5); + for (int dy = -1; dy <= 2; dy++) + { + float gy = max(0.0, abs(local.y - (float)dy) - 0.5); + for (int dx = -1; dx <= 2; dx++) + { + float gx = max(0.0, abs(local.x - (float)dx) - 0.5); + if (gx * gx + gy * gy + gz * gz >= 2.25) + continue; + + float3 gridStep = float3(dx, dy, dz); + int3 hash = cultmath_pcg3d(int3(cell + gridStep)); + float3 jitter = float3( + ((uint)hash.x >> 8) * (1.0 / 16777216.0) - 0.5, + ((uint)hash.y >> 8) * (1.0 / 16777216.0) - 0.5, + ((uint)hash.z >> 8) * (1.0 / 16777216.0) - 0.5); + float3 v = local - gridStep - jitter; + float falloff = exp(-2.0 * dot(v, v)); + float w = max(0.0, falloff - 0.01111); + float3 dw = w > 0.0 ? (-4.0 * falloff) * v : float3(0.0, 0.0, 0.0); + float phase = dot(v, side) + phaseOffset; + float c = cos(phase); + float s = sin(phase); + + sumW += w; + sumCos += w * c; + sumSin += w * s; + gradW += dw; + gradCos += dw * c - (w * s) * side; + gradSin += dw * s + (w * c) * side; + } + } + } + + float invW = 1.0 / sumW; + float rawCos = sumCos * invW; + float rawSin = sumSin * invW; + float3 gradRawCos = (gradCos - rawCos * gradW) * invW; + float3 gradRawSin = (gradSin - rawSin * gradW) * invW; + float rawLength = sqrt(rawCos * rawCos + rawSin * rawSin); + float floorLength = 1.0 - normalization; + + CultPhasor result; + if (rawLength > floorLength) + { + float outCos = rawCos / rawLength; + float outSin = rawSin / rawLength; + float3 gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength; + result.cos = float4((gradRawCos - outCos * gradLength) / rawLength, outCos); + result.sin = float4((gradRawSin - outSin * gradLength) / rawLength, outSin); + } + else + { + result.cos = float4(gradRawCos / floorLength, rawCos / floorLength); + result.sin = float4(gradRawSin / floorLength, rawSin / floorLength); + } + + return result; +} + +#endif diff --git a/packages/cultmath/shaders/CultMath.Phacelle.hlsl.meta b/packages/cultmath/shaders/CultMath.Phacelle.hlsl.meta new file mode 100644 index 00000000..69308465 --- /dev/null +++ b/packages/cultmath/shaders/CultMath.Phacelle.hlsl.meta @@ -0,0 +1,7 @@ +fileFormatVersion: 2 +guid: e07a3b95c1d84f2aa6b83d4c5f19a0d7 +ShaderIncludeImporter: + externalObjects: {} + userData: + assetBundleName: + assetBundleVariant: diff --git a/packages/cultmath/shaders/CultMath.hlsl b/packages/cultmath/shaders/CultMath.hlsl index d0dd22ae..af2d3741 100644 --- a/packages/cultmath/shaders/CultMath.hlsl +++ b/packages/cultmath/shaders/CultMath.hlsl @@ -503,91 +503,7 @@ CultCellular cultmath_cellular(float3 p) return result; } -// Invariant 8's struct return shape for cultmath_phacelle, compared field by field and bit for bit -// against CultMath.CultPhasor, same field order. -struct CultPhasor -{ - float4 cos; - float4 sin; -}; - -// Johansen's Phacelle noise on 3D cells. See math.phacelle's comment (math.cs) for the stripe wave -// vector, the weight, the normalization, the exact gradient and why the prune is exact. Mirrors -// the C# body operation for operation; the C#-only internal overload that switches the prune off -// has no shader twin. -CultPhasor cultmath_phacelle(float3 p, float3 side, float offset, float normalization) -{ - float3 cell = floor(p); - float3 local = p - cell; - float phaseOffset = offset * CULTMATH_TAU; - float sumW = 0.0; - float sumCos = 0.0; - float sumSin = 0.0; - float3 gradW = float3(0.0, 0.0, 0.0); - float3 gradCos = float3(0.0, 0.0, 0.0); - float3 gradSin = float3(0.0, 0.0, 0.0); - - for (int dz = -1; dz <= 2; dz++) - { - float gz = max(0.0, abs(local.z - (float)dz) - 0.5); - for (int dy = -1; dy <= 2; dy++) - { - float gy = max(0.0, abs(local.y - (float)dy) - 0.5); - for (int dx = -1; dx <= 2; dx++) - { - float gx = max(0.0, abs(local.x - (float)dx) - 0.5); - if (gx * gx + gy * gy + gz * gz >= 2.25) - continue; - - float3 gridStep = float3(dx, dy, dz); - int3 hash = cultmath_pcg3d(int3(cell + gridStep)); - float3 jitter = float3( - ((uint)hash.x >> 8) * (1.0 / 16777216.0) - 0.5, - ((uint)hash.y >> 8) * (1.0 / 16777216.0) - 0.5, - ((uint)hash.z >> 8) * (1.0 / 16777216.0) - 0.5); - float3 v = local - gridStep - jitter; - float falloff = exp(-2.0 * dot(v, v)); - float w = max(0.0, falloff - 0.01111); - float3 dw = w > 0.0 ? (-4.0 * falloff) * v : float3(0.0, 0.0, 0.0); - float phase = dot(v, side) + phaseOffset; - float c = cos(phase); - float s = sin(phase); - - sumW += w; - sumCos += w * c; - sumSin += w * s; - gradW += dw; - gradCos += dw * c - (w * s) * side; - gradSin += dw * s + (w * c) * side; - } - } - } - - float invW = 1.0 / sumW; - float rawCos = sumCos * invW; - float rawSin = sumSin * invW; - float3 gradRawCos = (gradCos - rawCos * gradW) * invW; - float3 gradRawSin = (gradSin - rawSin * gradW) * invW; - float rawLength = sqrt(rawCos * rawCos + rawSin * rawSin); - float floorLength = 1.0 - normalization; - - CultPhasor result; - if (rawLength > floorLength) - { - float outCos = rawCos / rawLength; - float outSin = rawSin / rawLength; - float3 gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength; - result.cos = float4((gradRawCos - outCos * gradLength) / rawLength, outCos); - result.sin = float4((gradRawSin - outSin * gradLength) / rawLength, outSin); - } - else - { - result.cos = float4(gradRawCos / floorLength, rawCos / floorLength); - result.sin = float4(gradRawSin / floorLength, rawSin / floorLength); - } - - return result; -} +#include "CultMath.Phacelle.hlsl" float cultmath_value_noise(float2 position) { diff --git a/packages/cultmath/src/CultMath/CultMath.csproj b/packages/cultmath/src/CultMath/CultMath.csproj index 7c253027..b4a6ce8b 100644 --- a/packages/cultmath/src/CultMath/CultMath.csproj +++ b/packages/cultmath/src/CultMath/CultMath.csproj @@ -8,6 +8,8 @@ CultMath CultMath GameCult.Math + + MIT AND MPL-2.0 0.2.4 Clean-room C# math shaped for HLSL-style renderer code. GameCult @@ -21,6 +23,8 @@ + + diff --git a/packages/cultmath/src/CultMath/CultPhasor.cs b/packages/cultmath/src/CultMath/CultPhasor.cs index 0a28faa8..fbc1f9d7 100644 --- a/packages/cultmath/src/CultMath/CultPhasor.cs +++ b/packages/cultmath/src/CultMath/CultPhasor.cs @@ -1,3 +1,9 @@ +// This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of +// the MPL was not distributed with this file, You can obtain one at https://mozilla.org/MPL/2.0/. +// +// The result type of math.phacelle (math.Phacelle.cs), Rune Skovbo Johansen's Phacelle noise +// (Copyright (c) 2025, MPL-2.0; see THIRD-PARTY-NOTICES.md). + namespace CultMath; /// diff --git a/packages/cultmath/src/CultMath/math.Phacelle.cs b/packages/cultmath/src/CultMath/math.Phacelle.cs new file mode 100644 index 00000000..25689737 --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Phacelle.cs @@ -0,0 +1,127 @@ +// This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of +// the MPL was not distributed with this file, You can obtain one at https://mozilla.org/MPL/2.0/. +// +// Phacelle noise. Copyright (c) 2025 Rune Skovbo Johansen (Shadertoy t3dyWl; MPL-2.0), read from the +// C# port PhacelleNoise in lpmitchell/AdvancedTerrainErosion, converted to burstable C# by Luke +// Mitchell (2026). Modified by GameCult: 3D cells, a caller-supplied stripe wave vector, an exact +// gradient and exact pruning. This file and CultPhasor.cs are the MPL-2.0 files of CultMath; the +// rest of the package is MIT (see THIRD-PARTY-NOTICES.md). + +namespace CultMath; + +public static partial class math +{ + // Johansen's Phacelle noise (phase + cell), generalised to 3D cells with a caller-supplied + // stripe wave vector. Each of the 4x4x4 cells around p carries one jittered feature point + // (jitter in [-0.5, 0.5) per axis, from pcg3d over the integer cell coordinate, the same + // integer-hashing rule as cellular) and one cosine/sine wave that runs along `side` (the wave + // vector: its length is 2*pi times the stripe frequency, and the caller chooses which + // perpendicular it is; on a unit sphere that is cross(p_hat, flow)). The waves are blended with + // the bell weight max(0, exp(-2 d^2) - 0.01111) (exactly 0 at distance 1.5, the nearest a + // feature outside the 4x4x4 block can be in 2D, which is what keeps cell-grid lines out), and + // the blended phasor is normalized: its length is stretched to 1 wherever the raw length is at + // least 1 - normalization, and scaled by 1 / (1 - normalization) below that. `offset` is the + // phase in cycles (1 is a full turn). Value and normalization follow Johansen's + // PhacelleNoise as ported in lpmitchell/AdvancedTerrainErosion. + // + // `normalization` is meant to lie in [0, 1]. At 0 the output is the raw blend (a phasor of + // length at most 1); at 1 the floor is 0 and every output is stretched to length 1. Values of + // 1 or more saturate at that (the floor is never above the raw length), and a negative value + // scales every output down by 1 / (1 - normalization) < 1. The normalized gradient scales as + // 1 / |I| where I is the blended phasor (below), so it grows wherever the blend nearly + // cancels; a normalization near 1 lets that happen at every point. + // + // The gradient is exact: it carries the weight derivative (grad w = -4 v exp(-2 |v|^2) where + // w > 0) and the normalization, not upstream's convention grad(cos) ~ -sin * side, which + // ignores the weight gradient (Asura probe: 16.8 deg off at the median, 97.8 at the 95th + // percentile; the exact form is 0.00 and 0.03). With P = sum(w * (cos, sin)) and W = sum(w), + // I = P / W has grad I = (grad P - I * grad W) / W; above the threshold the output is I / |I| + // with gradient (grad I - out * grad|I|) / |I|, below it grad I / (1 - normalization). The + // value is continuous at both seams below; the gradient is not. + // * Where |I| crosses the threshold, the gradient jumps between the two branches. + // * Where a cell's weight leaves its support (d^2 = ln(1 / 0.01111) / 2 = 2.24995, distance + // 1.49998), grad w jumps from -4 * 0.01111 * v to 0: a per-cell jump of length + // 0.0444 * |v|, 0.0667 at |v| = 1.5. Measured over random points, the resulting jump in + // the output gradient is a median 1.5% of its length and at most about 27%. + // + // W is the sum of 64 weights and is positive for every input pcg3d produces (the smallest + // W found by hill-climbing the hash-driven cell contents is 0.056). It is 0, and the output + // NaN, only if all 64 weights vanish, which takes adversarial jitter no hash yields. + // + // Pruning is exact: per axis, |local - g| - 0.5 (floored at 0) lower-bounds |v| for any jitter, + // so a cell whose summed square is at least 2.25 lies outside the support (2.24995 < 2.25) and + // has weight exactly 0: it contributes nothing to any sum. Cells whose bound sits just under + // 2.25 can still have a nonzero weight, so the cut cannot be lowered. On average the prune + // keeps about 45 of the 64 cells (a saving of about 30%); only about 14 have a nonzero weight. + // Only the internal overload can switch the prune off; PhacelleTests uses it to pin bit equality. + public static CultPhasor phacelle(float3 p, float3 side, float offset, float normalization) => + phacelle(p, side, offset, normalization, true); + + internal static CultPhasor phacelle(float3 p, float3 side, float offset, float normalization, bool prune) + { + var cell = floor(p); + var local = p - cell; + var phaseOffset = offset * TAU; + var sumW = 0.0f; + var sumCos = 0.0f; + var sumSin = 0.0f; + var gradW = new float3(0.0f, 0.0f, 0.0f); + var gradCos = new float3(0.0f, 0.0f, 0.0f); + var gradSin = new float3(0.0f, 0.0f, 0.0f); + + for (var dz = -1; dz <= 2; dz++) + { + var gz = max(0.0f, abs(local.z - (float)dz) - 0.5f); + for (var dy = -1; dy <= 2; dy++) + { + var gy = max(0.0f, abs(local.y - (float)dy) - 0.5f); + for (var dx = -1; dx <= 2; dx++) + { + var gx = max(0.0f, abs(local.x - (float)dx) - 0.5f); + if (prune && gx * gx + gy * gy + gz * gz >= 2.25f) + continue; + + var gridStep = new float3(dx, dy, dz); + var hash = pcg3d(int3(cell + gridStep)); + var jitter = new float3(cellular_unit(hash.x) - 0.5f, cellular_unit(hash.y) - 0.5f, cellular_unit(hash.z) - 0.5f); + var v = local - gridStep - jitter; + var falloff = exp(-2.0f * dot(v, v)); + var w = max(0.0f, falloff - 0.01111f); + var dw = w > 0.0f ? (-4.0f * falloff) * v : new float3(0.0f, 0.0f, 0.0f); + var phase = dot(v, side) + phaseOffset; + var c = cos(phase); + var s = sin(phase); + + sumW += w; + sumCos += w * c; + sumSin += w * s; + gradW += dw; + gradCos += dw * c - (w * s) * side; + gradSin += dw * s + (w * c) * side; + } + } + } + + var invW = 1.0f / sumW; + var rawCos = sumCos * invW; + var rawSin = sumSin * invW; + var gradRawCos = (gradCos - rawCos * gradW) * invW; + var gradRawSin = (gradSin - rawSin * gradW) * invW; + var rawLength = sqrt(rawCos * rawCos + rawSin * rawSin); + var floorLength = 1.0f - normalization; + + if (rawLength > floorLength) + { + var outCos = rawCos / rawLength; + var outSin = rawSin / rawLength; + var gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength; + return new CultPhasor( + new float4((gradRawCos - outCos * gradLength) / rawLength, outCos), + new float4((gradRawSin - outSin * gradLength) / rawLength, outSin)); + } + + return new CultPhasor( + new float4(gradRawCos / floorLength, rawCos / floorLength), + new float4(gradRawSin / floorLength, rawSin / floorLength)); + } +} diff --git a/packages/cultmath/src/CultMath/math.Phacelle.cs.meta b/packages/cultmath/src/CultMath/math.Phacelle.cs.meta new file mode 100644 index 00000000..898abc4c --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Phacelle.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 2c9d4e6b1a8f47d3b05e7a3c9d1f8e26 diff --git a/packages/cultmath/src/CultMath/math.cs b/packages/cultmath/src/CultMath/math.cs index 430325ea..23916b6b 100644 --- a/packages/cultmath/src/CultMath/math.cs +++ b/packages/cultmath/src/CultMath/math.cs @@ -968,104 +968,6 @@ public static partial class math // dxc and C#, bit-exact on both sides. internal static float cellular_unit(int bits) => ((uint)bits >> 8) * (1.0f / 16777216.0f); - // Johansen's Phacelle noise (phase + cell), generalised to 3D cells with a caller-supplied - // stripe wave vector. Each of the 4x4x4 cells around p carries one jittered feature point - // (jitter in [-0.5, 0.5) per axis, from pcg3d over the integer cell coordinate, the same - // integer-hashing rule as cellular) and one cosine/sine wave that runs along `side` (the wave - // vector: its length is 2*pi times the stripe frequency, and the caller chooses which - // perpendicular it is; on a unit sphere that is cross(p_hat, flow)). The waves are blended with - // the bell weight max(0, exp(-2 d^2) - 0.01111) (exactly 0 at distance 1.5, the nearest a - // feature outside the 4x4x4 block can be in 2D, which is what keeps cell-grid lines out), and - // the blended phasor is normalized: its length is stretched to 1 wherever the raw length is at - // least 1 - normalization, and scaled by 1 / (1 - normalization) below that. `offset` is the - // phase in cycles (1 is a full turn). Value and normalization follow Johansen's - // PhacelleNoise as ported in lpmitchell/AdvancedTerrainErosion. - // - // The gradient is exact: it carries the weight derivative (grad w = -4 v exp(-2 |v|^2) where - // w > 0) and the normalization, not upstream's convention grad(cos) ~ -sin * side, which - // ignores the weight gradient (Asura probe: 16.8 deg off at the median, 97.8 at the 95th - // percentile; the exact form is 0.00 and 0.03). With P = sum(w * (cos, sin)) and W = sum(w), - // I = P / W has grad I = (grad P - I * grad W) / W; above the threshold the output is I / |I| - // with gradient (grad I - out * grad|I|) / |I|, below it grad I / (1 - normalization). The - // gradient is discontinuous where |I| crosses the threshold, and where a cell's weight leaves - // its support (distance 1.4998): there grad w jumps from -4 * 0.01111 * v (length about 0.067) to 0. The value is - // continuous at both. - // - // Pruning is exact: per axis, |local - g| - 0.5 (floored at 0) lower-bounds |v| for any jitter, - // so a cell whose summed square is at least 2.25 has weight exactly 0 (exp(-4.5) = 0.011109 < - // 0.01111) and contributes nothing to any sum. Only the internal overload can switch the prune - // off; PhacelleTests uses it to pin bit equality. About 14 of 64 cells survive on average. - public static CultPhasor phacelle(float3 p, float3 side, float offset, float normalization) => - phacelle(p, side, offset, normalization, true); - - internal static CultPhasor phacelle(float3 p, float3 side, float offset, float normalization, bool prune) - { - var cell = floor(p); - var local = p - cell; - var phaseOffset = offset * TAU; - var sumW = 0.0f; - var sumCos = 0.0f; - var sumSin = 0.0f; - var gradW = new float3(0.0f, 0.0f, 0.0f); - var gradCos = new float3(0.0f, 0.0f, 0.0f); - var gradSin = new float3(0.0f, 0.0f, 0.0f); - - for (var dz = -1; dz <= 2; dz++) - { - var gz = max(0.0f, abs(local.z - (float)dz) - 0.5f); - for (var dy = -1; dy <= 2; dy++) - { - var gy = max(0.0f, abs(local.y - (float)dy) - 0.5f); - for (var dx = -1; dx <= 2; dx++) - { - var gx = max(0.0f, abs(local.x - (float)dx) - 0.5f); - if (prune && gx * gx + gy * gy + gz * gz >= 2.25f) - continue; - - var gridStep = new float3(dx, dy, dz); - var hash = pcg3d(int3(cell + gridStep)); - var jitter = new float3(cellular_unit(hash.x) - 0.5f, cellular_unit(hash.y) - 0.5f, cellular_unit(hash.z) - 0.5f); - var v = local - gridStep - jitter; - var falloff = exp(-2.0f * dot(v, v)); - var w = max(0.0f, falloff - 0.01111f); - var dw = w > 0.0f ? (-4.0f * falloff) * v : new float3(0.0f, 0.0f, 0.0f); - var phase = dot(v, side) + phaseOffset; - var c = cos(phase); - var s = sin(phase); - - sumW += w; - sumCos += w * c; - sumSin += w * s; - gradW += dw; - gradCos += dw * c - (w * s) * side; - gradSin += dw * s + (w * c) * side; - } - } - } - - var invW = 1.0f / sumW; - var rawCos = sumCos * invW; - var rawSin = sumSin * invW; - var gradRawCos = (gradCos - rawCos * gradW) * invW; - var gradRawSin = (gradSin - rawSin * gradW) * invW; - var rawLength = sqrt(rawCos * rawCos + rawSin * rawSin); - var floorLength = 1.0f - normalization; - - if (rawLength > floorLength) - { - var outCos = rawCos / rawLength; - var outSin = rawSin / rawLength; - var gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength; - return new CultPhasor( - new float4((gradRawCos - outCos * gradLength) / rawLength, outCos), - new float4((gradRawSin - outSin * gradLength) / rawLength, outSin)); - } - - return new CultPhasor( - new float4(gradRawCos / floorLength, rawCos / floorLength), - new float4(gradRawSin / floorLength, rawSin / floorLength)); - } - public static float value_noise(float2 position) { var cell = floor(position); diff --git a/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs index 2b695fcb..6b69a908 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs @@ -7,7 +7,7 @@ public sealed class HlslMirrorTests [Fact] public void HlslMirrorDoesNotPublishGeometryKernels() { - var source = File.ReadAllText(Path.Combine(GetCultMathRoot(), "shaders", "CultMath.hlsl")); + var source = HlslSourceCompatibilityTests.ReadShaderSource(GetCultMathRoot()); Assert.DoesNotContain("AdvancedErosionFilter.hlsl", source); Assert.DoesNotContain("Planetary.hlsl", source); Assert.DoesNotContain("spherical_erosion", source); @@ -17,7 +17,7 @@ public sealed class HlslMirrorTests [Fact] public void HlslMirrorPublishesCultMathPrimitives() { - var include = File.ReadAllText(Path.Combine(GetCultMathRoot(), "shaders", "CultMath.hlsl")); + var include = HlslSourceCompatibilityTests.ReadShaderSource(GetCultMathRoot()); var requiredSymbols = new[] { "cultmath_radians", @@ -62,6 +62,17 @@ public sealed class HlslMirrorTests Assert.DoesNotContain("cultmath_normalize", include); } + [Theory] + [InlineData("CultMath.hlsl")] + [InlineData("CultMath.Phacelle.hlsl")] + public void UnityPackageShaderIsIdenticalToTheCanonicalShader(string name) + { + var root = GetCultMathRoot(); + Assert.Equal( + File.ReadAllText(Path.Combine(root, "shaders", name)).ReplaceLineEndings("\n"), + File.ReadAllText(Path.Combine(root, "unity", "org.gamecult.cultmath", "Shaders", name)).ReplaceLineEndings("\n")); + } + private static string GetCultMathRoot() { var directory = new DirectoryInfo(AppContext.BaseDirectory); diff --git a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs index 2cf737df..2cb8ae49 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs @@ -163,6 +163,33 @@ public sealed class HlslSourceCompatibilityTests } } + /// + /// The dense sweep almost never meets a cell whose box bound is just under the 2.25 cut and whose + /// weight is still nonzero: a mirror cut of 2.245, 2.24 or 2.2 survives it. These points sit next + /// to such cells (), so the mirror must match C# bit for + /// bit, and must match the unpruned C# sum (the shader has no prune switch, so that is the + /// independent statement that the cut is exact). + /// + [Fact] + public void PhacelleMirrorMatchesCSharpBitForBitNextToCellsJustInsideTheCut() + { + var (assembly, errors) = CompileShaderMirror(); + Assert.True(assembly is not null, string.Join(Environment.NewLine, errors)); + var shaderType = assembly!.GetType("CultMathHlsl.HlslShader")!; + var shader = Activator.CreateInstance(shaderType); + var mirror = ShaderFunction(shaderType, "cultmath_phacelle", typeof(float3), typeof(float3), typeof(float), typeof(float))!; + + foreach (var p in PhacelleTests.PruneCutPoints) + { + Assert.True(PhacelleTests.HasCellJustInsideTheCut(p), $"no cell just inside the cut at {p}"); + var actual = Values(mirror.Invoke(shader, new object[] { p, PhacelleTests.PruneCutSide, PhacelleTests.PruneCutOffset, PhacelleTests.PruneCutNormalization })!); + var pruned = Values(phacelle(p, PhacelleTests.PruneCutSide, PhacelleTests.PruneCutOffset, PhacelleTests.PruneCutNormalization)); + var unpruned = Values(phacelle(p, PhacelleTests.PruneCutSide, PhacelleTests.PruneCutOffset, PhacelleTests.PruneCutNormalization, false)); + Assert.True(BitwiseEqual(pruned, actual), $"mirror differs from C# at p = {p}"); + Assert.True(BitwiseEqual(unpruned, actual), $"mirror differs from the unpruned C# sum at p = {p}"); + } + } + /// /// A hand-built pair of structs standing in for a mirror struct return and its C# counterpart /// (invariant 8's one struct return shape, e.g. ): same field shape, @@ -283,6 +310,18 @@ public sealed class HlslSourceCompatibilityTests ? (float.IsNaN(a) && float.IsNaN(b)) || BitConverter.SingleToInt32Bits(a) == BitConverter.SingleToInt32Bits(b) : p.First.Value.Equals(p.Second.Value))); + /// + /// The text of CultMath.hlsl with each #include "name" replaced in place by the + /// included file from the same directory, as the shader compiler would resolve it. Only the + /// package's own includes exist, so a missing file is an error rather than a skip. + /// + internal static string ReadShaderSource(string cultMathRoot) + { + var shaders = Path.Combine(cultMathRoot, "shaders"); + return Regex.Replace(File.ReadAllText(Path.Combine(shaders, "CultMath.hlsl")), @"(?m)^#include ""([^""]+)""[ \t]*\r?$", + match => File.ReadAllText(Path.Combine(shaders, match.Groups[1].Value))); + } + /// The documented HLSL-to-C# transformations, and nothing else. internal static string TransformHlslToCSharp(string hlsl) { @@ -311,7 +350,7 @@ public sealed class HlslSourceCompatibilityTests private static (Assembly? Assembly, IReadOnlyList Errors) CompileShaderMirror() { - var hlsl = File.ReadAllText(Path.Combine(FindCultMathRoot(), "shaders", "CultMath.hlsl")); + var hlsl = ReadShaderSource(FindCultMathRoot()); var tree = CSharpSyntaxTree.ParseText(TransformHlslToCSharp(hlsl)); var references = ((string)AppContext.GetData("TRUSTED_PLATFORM_ASSEMBLIES")!) .Split(Path.PathSeparator) diff --git a/packages/cultmath/tests/CultMath.Tests/PhacelleTests.cs b/packages/cultmath/tests/CultMath.Tests/PhacelleTests.cs index 15d1a89a..b64f70b1 100644 --- a/packages/cultmath/tests/CultMath.Tests/PhacelleTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/PhacelleTests.cs @@ -144,7 +144,7 @@ public sealed class PhacelleTests { // Above the threshold the output is I / |I|; below it, I / (1 - normalization). The gradient // is discontinuous where |I| crosses the threshold, and at the edge of each cell's weight support - // (distance 1.4998, where max(0, ...) starts: a jump of about 4 d 0.01111 in grad w), so a band + // (distance 1.49998, where max(0, ...) starts: a jump of 4 d 0.01111 in grad w), so a band // around each is excluded (by the // reference's raw length, since the output alone cannot say how close it is). const float normalization = 0.5f; @@ -265,6 +265,71 @@ public sealed class PhacelleTests } } + // Points next to a cell whose jitter sits at the far corner of its cube: the cell's box bound is + // just under the 2.25 cut (2.2451 to 2.25) while its real distance is still under the support + // radius (d^2 < 2.24995), so its weight is small but nonzero and a prune that cuts even slightly + // low (2.245, 2.24, 2.2) drops a contribution and changes output bits. Random points meet such + // a cell about once in 10^4 evaluations, which is why the random sweep above let those cuts + // through. Each was found once by an offline search (cells K in [-12, 12]^3 whose jitter on one + // axis lies within 0.0008 of the cube face; p is the feature moved 1.4996 out along that axis, + // so |p - feature| = 1.4996 and the bound is about 1.4988). The witness check below re-derives + // that property from the hash, so a change to pcg3d fails loudly instead of degenerating. + // Shared with the HLSL mirror test. + internal static readonly float3[] PruneCutPoints = + { + new(1.77850008f, -2.16020727f, -0.9991166f), + new(1.09346282f, -2.37882328f, -3.00082445f), + new(-6.39144945f, -3.99953747f, -0.610023975f), + new(-0.99940908f, 4.16302729f, 4.07277536f), + new(3.99912024f, 7.21890545f, 1.09053993f), + new(5.43959713f, 1.99937367f, -6.14929914f), + }; + + internal static readonly float3 PruneCutSide = new(2.1f, -3.7f, 1.3f); + internal const float PruneCutOffset = 0.3f; + internal const float PruneCutNormalization = 0.5f; + + // True when some cell of the 4x4x4 window has a box bound in [2.2451, 2.25) and a positive weight. + internal static bool HasCellJustInsideTheCut(float3 p) + { + double fx = Math.Floor(p.x), fy = Math.Floor(p.y), fz = Math.Floor(p.z); + double lx = p.x - fx, ly = p.y - fy, lz = p.z - fz; + for (var dz = -1; dz <= 2; dz++) + for (var dy = -1; dy <= 2; dy++) + for (var dx = -1; dx <= 2; dx++) + { + var gx = Math.Max(0.0, Math.Abs(lx - dx) - 0.5); + var gy = Math.Max(0.0, Math.Abs(ly - dy) - 0.5); + var gz = Math.Max(0.0, Math.Abs(lz - dz) - 0.5); + var bound = gx * gx + gy * gy + gz * gz; + if (bound < 2.2451 || bound >= 2.25) + continue; + + var hash = math.pcg3d(new int3((int)fx + dx, (int)fy + dy, (int)fz + dz)); + double vx = lx - dx - (math.cellular_unit(hash.x) - 0.5); + double vy = ly - dy - (math.cellular_unit(hash.y) - 0.5); + double vz = lz - dz - (math.cellular_unit(hash.z) - 0.5); + if (Math.Exp(-2.0 * (vx * vx + vy * vy + vz * vz)) > 0.01111) + return true; + } + + return false; + } + + [Fact] + public void PrunedOutputIsBitEqualToTheUnprunedSumNextToCellsJustInsideTheCut() + { + foreach (var p in PruneCutPoints) + { + Assert.True(HasCellJustInsideTheCut(p), $"no cell just inside the cut at {p}: the fixture no longer exercises the prune"); + var pruned = math.phacelle(p, PruneCutSide, PruneCutOffset, PruneCutNormalization); + var unpruned = math.phacelle(p, PruneCutSide, PruneCutOffset, PruneCutNormalization, false); + + AssertBitEqual(unpruned.cos, pruned.cos); + AssertBitEqual(unpruned.sin, pruned.sin); + } + } + private static void AssertBitEqual(float4 expected, float4 actual) { Assert.Equal(BitConverter.SingleToInt32Bits(expected.x), BitConverter.SingleToInt32Bits(actual.x)); diff --git a/packages/cultmath/unity/org.gamecult.cultmath/LICENSE.md b/packages/cultmath/unity/org.gamecult.cultmath/LICENSE.md index 5e50f2d6..7e18d22d 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/LICENSE.md +++ b/packages/cultmath/unity/org.gamecult.cultmath/LICENSE.md @@ -4,7 +4,7 @@ CultMath is licensed under the Mozilla Public License 2.0 (`MPL-2.0`). The complete license text is available from the Mozilla Foundation at . Source and modification history for the -precompiled assembly in this package are available in the GameCult/CultMath -repository at the corresponding package tag. +precompiled assembly in this package are available in the GameCult/CultLib +repository (`packages/cultmath`) at the corresponding package tag. Third-party code retains its original license. See `THIRD-PARTY-NOTICES.md`. diff --git a/packages/cultmath/unity/org.gamecult.cultmath/README.md b/packages/cultmath/unity/org.gamecult.cultmath/README.md index 20e1840c..4981c54c 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/README.md +++ b/packages/cultmath/unity/org.gamecult.cultmath/README.md @@ -2,7 +2,7 @@ This is CultMath's Unity 2021.3-compatible Git-UPM surface. It contains the precompiled `CultMath.dll`, portable symbols, the numeric -`Shaders/CultMath.hlsl` mirror, and the `CultMath.UnityBridge` assembly. +`Shaders/CultMath.hlsl` mirror (with its `CultMath.Phacelle.hlsl` include), and the `CultMath.UnityBridge` assembly. `CultMath.UnityBridge` is the only C# source Unity compiles here. It holds the UnityEngine conversions the engine-free core cannot declare: `ToCultMath()` on diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Phacelle.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Phacelle.hlsl new file mode 100644 index 00000000..1db0a2fc --- /dev/null +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Phacelle.hlsl @@ -0,0 +1,99 @@ +#ifndef CULTMATH_PHACELLE_HLSL +#define CULTMATH_PHACELLE_HLSL + +// This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of +// the MPL was not distributed with this file, You can obtain one at https://mozilla.org/MPL/2.0/. +// +// Phacelle noise. Copyright (c) 2025 Rune Skovbo Johansen (Shadertoy t3dyWl; MPL-2.0), read from the +// C# port PhacelleNoise in lpmitchell/AdvancedTerrainErosion, converted to burstable C# by Luke +// Mitchell (2026). Modified by GameCult: 3D cells, a caller-supplied stripe wave vector, an exact +// gradient and exact pruning. This file is an MPL-2.0 file; CultMath.hlsl, which includes it after +// its own definitions of CULTMATH_TAU and cultmath_pcg3d, is MIT (see THIRD-PARTY-NOTICES.md). + +// Invariant 8's struct return shape for cultmath_phacelle, compared field by field and bit for bit +// against CultMath.CultPhasor, same field order. +struct CultPhasor +{ + float4 cos; + float4 sin; +}; + +// Johansen's Phacelle noise on 3D cells. See math.phacelle's comment (math.Phacelle.cs) for the stripe wave +// vector, the weight, the normalization, the exact gradient and why the prune is exact. Mirrors +// the C# body operation for operation; the C#-only internal overload that switches the prune off +// has no shader twin. +CultPhasor cultmath_phacelle(float3 p, float3 side, float offset, float normalization) +{ + float3 cell = floor(p); + float3 local = p - cell; + float phaseOffset = offset * CULTMATH_TAU; + float sumW = 0.0; + float sumCos = 0.0; + float sumSin = 0.0; + float3 gradW = float3(0.0, 0.0, 0.0); + float3 gradCos = float3(0.0, 0.0, 0.0); + float3 gradSin = float3(0.0, 0.0, 0.0); + + for (int dz = -1; dz <= 2; dz++) + { + float gz = max(0.0, abs(local.z - (float)dz) - 0.5); + for (int dy = -1; dy <= 2; dy++) + { + float gy = max(0.0, abs(local.y - (float)dy) - 0.5); + for (int dx = -1; dx <= 2; dx++) + { + float gx = max(0.0, abs(local.x - (float)dx) - 0.5); + if (gx * gx + gy * gy + gz * gz >= 2.25) + continue; + + float3 gridStep = float3(dx, dy, dz); + int3 hash = cultmath_pcg3d(int3(cell + gridStep)); + float3 jitter = float3( + ((uint)hash.x >> 8) * (1.0 / 16777216.0) - 0.5, + ((uint)hash.y >> 8) * (1.0 / 16777216.0) - 0.5, + ((uint)hash.z >> 8) * (1.0 / 16777216.0) - 0.5); + float3 v = local - gridStep - jitter; + float falloff = exp(-2.0 * dot(v, v)); + float w = max(0.0, falloff - 0.01111); + float3 dw = w > 0.0 ? (-4.0 * falloff) * v : float3(0.0, 0.0, 0.0); + float phase = dot(v, side) + phaseOffset; + float c = cos(phase); + float s = sin(phase); + + sumW += w; + sumCos += w * c; + sumSin += w * s; + gradW += dw; + gradCos += dw * c - (w * s) * side; + gradSin += dw * s + (w * c) * side; + } + } + } + + float invW = 1.0 / sumW; + float rawCos = sumCos * invW; + float rawSin = sumSin * invW; + float3 gradRawCos = (gradCos - rawCos * gradW) * invW; + float3 gradRawSin = (gradSin - rawSin * gradW) * invW; + float rawLength = sqrt(rawCos * rawCos + rawSin * rawSin); + float floorLength = 1.0 - normalization; + + CultPhasor result; + if (rawLength > floorLength) + { + float outCos = rawCos / rawLength; + float outSin = rawSin / rawLength; + float3 gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength; + result.cos = float4((gradRawCos - outCos * gradLength) / rawLength, outCos); + result.sin = float4((gradRawSin - outSin * gradLength) / rawLength, outSin); + } + else + { + result.cos = float4(gradRawCos / floorLength, rawCos / floorLength); + result.sin = float4(gradRawSin / floorLength, rawSin / floorLength); + } + + return result; +} + +#endif diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Phacelle.hlsl.meta b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Phacelle.hlsl.meta new file mode 100644 index 00000000..9c074cb5 --- /dev/null +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Phacelle.hlsl.meta @@ -0,0 +1,7 @@ +fileFormatVersion: 2 +guid: 5b1c7e0a94d24f6c8a3e2d17f60b9c48 +ShaderIncludeImporter: + externalObjects: {} + userData: + assetBundleName: + assetBundleVariant: diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl index d0dd22ae..af2d3741 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl @@ -503,91 +503,7 @@ CultCellular cultmath_cellular(float3 p) return result; } -// Invariant 8's struct return shape for cultmath_phacelle, compared field by field and bit for bit -// against CultMath.CultPhasor, same field order. -struct CultPhasor -{ - float4 cos; - float4 sin; -}; - -// Johansen's Phacelle noise on 3D cells. See math.phacelle's comment (math.cs) for the stripe wave -// vector, the weight, the normalization, the exact gradient and why the prune is exact. Mirrors -// the C# body operation for operation; the C#-only internal overload that switches the prune off -// has no shader twin. -CultPhasor cultmath_phacelle(float3 p, float3 side, float offset, float normalization) -{ - float3 cell = floor(p); - float3 local = p - cell; - float phaseOffset = offset * CULTMATH_TAU; - float sumW = 0.0; - float sumCos = 0.0; - float sumSin = 0.0; - float3 gradW = float3(0.0, 0.0, 0.0); - float3 gradCos = float3(0.0, 0.0, 0.0); - float3 gradSin = float3(0.0, 0.0, 0.0); - - for (int dz = -1; dz <= 2; dz++) - { - float gz = max(0.0, abs(local.z - (float)dz) - 0.5); - for (int dy = -1; dy <= 2; dy++) - { - float gy = max(0.0, abs(local.y - (float)dy) - 0.5); - for (int dx = -1; dx <= 2; dx++) - { - float gx = max(0.0, abs(local.x - (float)dx) - 0.5); - if (gx * gx + gy * gy + gz * gz >= 2.25) - continue; - - float3 gridStep = float3(dx, dy, dz); - int3 hash = cultmath_pcg3d(int3(cell + gridStep)); - float3 jitter = float3( - ((uint)hash.x >> 8) * (1.0 / 16777216.0) - 0.5, - ((uint)hash.y >> 8) * (1.0 / 16777216.0) - 0.5, - ((uint)hash.z >> 8) * (1.0 / 16777216.0) - 0.5); - float3 v = local - gridStep - jitter; - float falloff = exp(-2.0 * dot(v, v)); - float w = max(0.0, falloff - 0.01111); - float3 dw = w > 0.0 ? (-4.0 * falloff) * v : float3(0.0, 0.0, 0.0); - float phase = dot(v, side) + phaseOffset; - float c = cos(phase); - float s = sin(phase); - - sumW += w; - sumCos += w * c; - sumSin += w * s; - gradW += dw; - gradCos += dw * c - (w * s) * side; - gradSin += dw * s + (w * c) * side; - } - } - } - - float invW = 1.0 / sumW; - float rawCos = sumCos * invW; - float rawSin = sumSin * invW; - float3 gradRawCos = (gradCos - rawCos * gradW) * invW; - float3 gradRawSin = (gradSin - rawSin * gradW) * invW; - float rawLength = sqrt(rawCos * rawCos + rawSin * rawSin); - float floorLength = 1.0 - normalization; - - CultPhasor result; - if (rawLength > floorLength) - { - float outCos = rawCos / rawLength; - float outSin = rawSin / rawLength; - float3 gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength; - result.cos = float4((gradRawCos - outCos * gradLength) / rawLength, outCos); - result.sin = float4((gradRawSin - outSin * gradLength) / rawLength, outSin); - } - else - { - result.cos = float4(gradRawCos / floorLength, rawCos / floorLength); - result.sin = float4(gradRawSin / floorLength, rawSin / floorLength); - } - - return result; -} +#include "CultMath.Phacelle.hlsl" float cultmath_value_noise(float2 position) { diff --git a/packages/cultmath/unity/org.gamecult.cultmath/THIRD-PARTY-NOTICES.md b/packages/cultmath/unity/org.gamecult.cultmath/THIRD-PARTY-NOTICES.md index 73f2461a..625d63a9 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/THIRD-PARTY-NOTICES.md +++ b/packages/cultmath/unity/org.gamecult.cultmath/THIRD-PARTY-NOTICES.md @@ -32,5 +32,10 @@ gradient and the exact pruning are new here. - Copyright (c) 2025 Rune Skovbo Johansen. - Upstream: -- License: MPL-2.0 (). CultMath's Unity package - is MPL-2.0, so the licence is compatible. +- License: MPL-2.0 (). MPL-2.0 is file-level, so it + covers exactly the files that carry the port and nothing else: + `src/CultMath/math.Phacelle.cs`, `src/CultMath/CultPhasor.cs`, + `shaders/CultMath.Phacelle.hlsl` and the Unity copy of that include. Each + carries the MPL-2.0 header and this attribution. The rest of CultMath is MIT, + so the `GameCult.Math` NuGet package declares `MIT AND MPL-2.0` and ships this + file. The Unity package as a whole is MPL-2.0.