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++)