From 1c913f5c6d978ffe16dcbc85eef2fb64073632c1 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 17 Oct 2026 06:44:07 +0200 Subject: [PATCH] cultmath: weighted_lattice_cells with exact ring-exit search, four public steps, range constants; power_border splits a private unchecked core Co-Authored-By: Claude Sonnet 5.5 --- .../cultmath/src/CultMath/math.Lattice.cs | 202 ++++++++++++++++++ .../src/CultMath/math.Lattice.cs.meta | 2 + packages/cultmath/src/CultMath/math.Sites.cs | 14 +- 3 files changed, 217 insertions(+), 1 deletion(-) create mode 100644 packages/cultmath/src/CultMath/math.Lattice.cs create mode 100644 packages/cultmath/src/CultMath/math.Lattice.cs.meta diff --git a/packages/cultmath/src/CultMath/math.Lattice.cs b/packages/cultmath/src/CultMath/math.Lattice.cs new file mode 100644 index 00000000..0dd0be56 --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Lattice.cs @@ -0,0 +1,202 @@ +namespace CultMath; + +public static partial class math +{ + // Weighted lattice cells: a jittered, drifting cellular field whose cells are weighted by how close the cell's + // centre lies to a border of the body diagram (the power diagram of the zone's body span). The public steps are + // the four per-cell rules a tape or a shader composes; weighted_lattice_cells is the loop over them. + + /// The greatest boost accepts: weights span 0.5 to 5, a ratio of 10:1. + public const float LatticeMaxBoost = 4.0f; + + /// The greatest stretch accepts. + public const float LatticeMaxStretch = 2.0f; + + /// The greatest smoothing accepts, in cell units. + public const float LatticeMaxSmoothing = 1.0f; + + /// + /// The outermost ring the search can reach anywhere in the accepted range (see the proof above + /// ). + /// + public const int LatticeMaxRing = 37; + + // A query beyond this many cells from the origin has no integer cell id in float precision. + private const float LatticeMaxCoordinate = 8388608.0f; + + /// + /// The site of a lattice cell: float3(site.xy, base weight). The site is cell + 0.5 + 0.5 sin(drift + /// time + 2 pi u) per axis, so it stays within one cell of the cell's corner and drifts at rate + /// . The base weight is 0.5 + 0.5 u, in [0.5, 1), from its own hash lane. + /// + public static float3 lattice_site(float2 cell, float time, float drift) + { + var h = pcg3d(new int3((int)cell.x, (int)cell.y, 0)); + var phase = new float2(cellular_unit(h.x), cellular_unit(h.y)); + var site = cell + 0.5f + 0.5f * new float2(sin(drift * time + TAU * phase.x), sin(drift * time + TAU * phase.y)); + return new float3(site, 0.5f + 0.5f * cellular_unit(h.z)); + } + + /// + /// The weight of a cell: base (1 + boost exp(-(border / falloff)^2)), where is + /// the distance from the cell's centre to the nearest body-diagram border. + /// + public static float lattice_weight(float @base, float border, float boost, float falloff) + { + var t = border / falloff; + return @base * (1.0f + boost * exp(-(t * t))); + } + + /// + /// The weighted distance from to , in cell units: + /// weight |(stretch a, c)| with a and c the offset's components along and across the unit + /// . Isotropic at stretch 1. + /// + public static float lattice_distance(float2 q, float2 site, float weight, float2 normal, float stretch) + { + var v = site - q; + var a = dot(v, normal); + var c = v.x * normal.y - v.y * normal.x; + return weight * length(new float2(stretch * a, c)); + } + + /// + /// One step of the smooth minimum fold: folded with distance , for a + /// nearest distance . A distance of at least f1 + 2 smoothing leaves + /// unchanged, and the result never exceeds f1 + smoothing. Start the fold at + /// f1 + smoothing and visit every cell once, in an order that does not depend on the search box. + /// + public static float lattice_smooth_step(float m, float d, float f1, float smoothing) + { + if (d >= f1 + 2.0f * smoothing) + return m; + var h = smoothstep(0.0f, 1.0f, 0.5f + 0.5f * (m - d) / smoothing); + var next = m + (d - m) * h - h * (1.0f - h) * smoothing / (1.0f + 3.0f * smoothing); + return min(next, f1 + smoothing); + } + + // Proof that the ring search finds every cell that can change an output (ruling lattice-ring-exit-sound). + // + // Let c0 = floor(q) and let a cell lie on Chebyshev ring k around c0 (an offset of k on some axis). Its site lies + // within [cell, cell + 1] per axis and q within [c0, c0 + 1), so on that axis |site - q| >= k - 1. The length of + // the offset is at least that, the stretch is at least 1 (it only lengthens), and the weight is at least + // 0.5 (base >= 0.5, boost >= 0). So every cell on ring k or beyond has d >= 0.5 (k - 1), and the same holds for + // every farther ring. + // + // A cell with d >= F2 changes neither F1 nor F2, and a cell with d >= F1 + 2 smoothing leaves the smooth minimum + // unchanged (lattice_smooth_step returns m before folding it). So the search stops before ring k once + // 0.5 (k - 1) >= max(F2, F1 + 2 smoothing), with F1 and F2 taken over the rings visited. Those values only fall + // as rings are added, so the running cutoff is at least the final one: stopping on it never stops early. Every + // unvisited cell then changes no output, for every weight, stretch and drift. + // + // The loop terminates inside the accepted range. Cell c0 gives a distance below sqrt(2) in cell units. The + // neighbour on the nearer side of the axis whose query fraction sits closest to a cell boundary (fraction m <= 1/2) + // lies within (1 + m, 1) per axis, below sqrt(3.25). Weights are below 1 + LatticeMaxBoost and stretch is at most + // LatticeMaxStretch, so F2 < 10 sqrt(3.25) = 18.03 and F1 + 2 smoothing < 10 sqrt(2) + 2 = 16.2. The exit holds at + // 0.5 (k - 1) >= 18.03, that is at k = 38, so rings 0 to 37 (LatticeMaxRing) are the most visited. + // + // The result does not depend on the box visited. The second pass folds F1, F2 and the winning id in absolute + // row-major order (strict <, so an exact tie keeps the earlier cell), and the smooth minimum starts at + // F1 + smoothing and is clamped to it, so a cell with d >= F1 + 2 smoothing, visited or not, changes nothing. + + /// + /// Weighted cellular noise over a jittered lattice: float4(F1, F2, smooth, id). F1 and F2 are the two + /// smallest weighted distances from in cell units times weight, smooth is the smooth minimum + /// (, in the same units), and id is a hash of the winning cell in [0, 1). + /// + /// + /// A cell's weight is of the exact power border distance () + /// of the cell's centre in the body span (float4(position.xyz, sigma), z = 0), so it is + /// continuous as the bodies move. Distances are stretched by along + /// (isotropic at 1). The search is exact in the accepted range: boost 0 to + /// , stretch 1 to , smoothing up to . + /// + /// An argument is not finite, a size is not above zero, a value is outside its range, the normal is zero while stretch exceeds 1, or the body span is refused as refuses it. + public static float4 weighted_lattice_cells( + float2 p, float time, ReadOnlySpan bodies, float cellSize, float drift, float boost, float falloff, + float smoothing, float2 stretchNormal, float stretch) + { + if (!float.IsFinite(p.x) || !float.IsFinite(p.y) || !float.IsFinite(cellSize) || !(cellSize > 0.0f) + || !(abs(p.x / cellSize) < LatticeMaxCoordinate) || !(abs(p.y / cellSize) < LatticeMaxCoordinate)) + throw new ArgumentException("p must be finite and within 2^23 cells of the origin, and cellSize finite and greater than zero.", nameof(p)); + if (!float.IsFinite(time)) + throw new ArgumentException("time must be finite.", nameof(time)); + if (!float.IsFinite(drift)) + throw new ArgumentException("drift must be finite.", nameof(drift)); + if (!float.IsFinite(boost) || boost < 0.0f || boost > LatticeMaxBoost) + throw new ArgumentException("boost must be within [0, LatticeMaxBoost].", nameof(boost)); + if (!float.IsFinite(falloff) || !(falloff > 0.0f)) + throw new ArgumentException("falloff must be finite and greater than zero.", nameof(falloff)); + if (!float.IsFinite(smoothing) || !(smoothing > 0.0f) || smoothing > LatticeMaxSmoothing) + throw new ArgumentException("smoothing must be within (0, LatticeMaxSmoothing].", nameof(smoothing)); + if (!float.IsFinite(stretch) || stretch < 1.0f || stretch > LatticeMaxStretch) + throw new ArgumentException("stretch must be within [1, LatticeMaxStretch].", nameof(stretch)); + if (!float.IsFinite(stretchNormal.x) || !float.IsFinite(stretchNormal.y)) + throw new ArgumentException("stretchNormal must be finite.", nameof(stretchNormal)); + var normalLength = length(stretchNormal); + if (normalLength > 0.0f) + stretchNormal /= normalLength; + else if (stretch > 1.0f) + throw new ArgumentException("stretchNormal must not be zero while stretch exceeds 1.", nameof(stretchNormal)); + else + stretchNormal = new float2(1.0f, 0.0f); + require_cell_span(new float3(p, 0.0f), bodies, positiveWeights: false); + + const int Width = 2 * LatticeMaxRing + 1; + Span cache = stackalloc float[Width * Width]; + var q = p / cellSize; + var origin = floor(q); + var cx0 = (int)origin.x; + var cz0 = (int)origin.y; + var f1 = 1.0e30f; + var f2 = 1.0e30f; + var reach = 0; + + for (var k = 0; ; k++) + { + if (k > 1 && 0.5f * (k - 1) >= max(f2, f1 + 2.0f * smoothing)) + break; + if (k > LatticeMaxRing) + throw new InvalidOperationException("The lattice ring search left its proven range."); + for (var dz = -k; dz <= k; dz++) + { + var step = dz == -k || dz == k || k == 0 ? 1 : 2 * k; + for (var dx = -k; dx <= k; dx += step) + { + var d = lattice_cell_distance(q, new float2(cx0 + dx, cz0 + dz), time, drift, bodies, cellSize, boost, falloff, stretchNormal, stretch); + cache[(dz + LatticeMaxRing) * Width + dx + LatticeMaxRing] = d; + if (d < f1) { f2 = f1; f1 = d; } + else if (d < f2) f2 = d; + } + } + reach = k; + } + + var m = f1 + smoothing; + var best = 1.0e30f; + var second = 1.0e30f; + var winDx = 0; + var winDz = 0; + for (var dz = -reach; dz <= reach; dz++) + for (var dx = -reach; dx <= reach; dx++) + { + var d = cache[(dz + LatticeMaxRing) * Width + dx + LatticeMaxRing]; + if (d < best) { second = best; best = d; winDx = dx; winDz = dz; } + else if (d < second) second = d; + m = lattice_smooth_step(m, d, f1, smoothing); + } + + var id = pcg4d(new int4(cx0 + winDx, cz0 + winDz, 0, 0)); + return new float4(best, second, m, cellular_unit(id.w)); + } + + private static float lattice_cell_distance( + float2 q, float2 cell, float time, float drift, ReadOnlySpan bodies, float cellSize, float boost, float falloff, + float2 normal, float stretch) + { + var site = lattice_site(cell, time, drift); + var centre = (cell + 0.5f) * cellSize; + var border = power_border_unchecked(new float3(centre, 0.0f), bodies).w; + return lattice_distance(q, site.xy, lattice_weight(site.z, border, boost, falloff), normal, stretch); + } +} diff --git a/packages/cultmath/src/CultMath/math.Lattice.cs.meta b/packages/cultmath/src/CultMath/math.Lattice.cs.meta new file mode 100644 index 00000000..a4f24ecd --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Lattice.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 0a182e056c8c4258881c347ead6d571d diff --git a/packages/cultmath/src/CultMath/math.Sites.cs b/packages/cultmath/src/CultMath/math.Sites.cs index feab825c..2509ec40 100644 --- a/packages/cultmath/src/CultMath/math.Sites.cs +++ b/packages/cultmath/src/CultMath/math.Sites.cs @@ -122,6 +122,11 @@ public static partial class math public static CultSiteCells power_cells(float3 p, ReadOnlySpan sites) { require_cell_span(p, sites, positiveWeights: false); + return power_cells_unchecked(p, sites); + } + + private static CultSiteCells power_cells_unchecked(float3 p, ReadOnlySpan sites) + { var cells = site_cells_empty(); for (var i = 0; i < sites.Length; i++) cells = site_cells_insert(cells, power_score(p, sites[i]), i); @@ -137,9 +142,16 @@ public static partial class math /// /// As . public static float4 power_border(float3 p, ReadOnlySpan sites) + { + require_cell_span(p, sites, positiveWeights: false); + return power_border_unchecked(p, sites); + } + + // The border over an already validated span, which the lattice calls once per cell. + private static float4 power_border_unchecked(float3 p, ReadOnlySpan sites) { var border = new float4(0.0f, 0.0f, 0.0f, 1.0e30f); - var owner = (int)power_cells(p, sites).nearestSite; + var owner = (int)power_cells_unchecked(p, sites).nearestSite; if (owner < 0) return border; for (var k = 0; k < sites.Length; k++)