From 7a656a7ec2b1df2bc2cad6ec506be94734ef489e Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Fri, 25 Sep 2026 18:06:29 +0200 Subject: [PATCH] Add surface nets extraction and face-weighted normals to GameCult.Geometry Implements CultGeometrySurfaceNets.Extract (welded CultGeometryQuadMesh, one quad per interior crossing grid edge, outward-oriented winding) and CultGeometryQuadNormals.FaceWeighted (unnormalized diagonal-cross-product area weighting, matching CultGeometryIsoSurface's sign convention and degenerate-normal rule). Fresh implementation per SURFACE-NETS-PROVENANCE.md; no code or tables copied from the cited literature. Co-Authored-By: Claude Sonnet 5 --- src/GameCult.Geometry/CultGeometryQuadMesh.cs | 88 +++++ .../CultGeometryQuadNormals.cs | 65 ++++ .../CultGeometrySurfaceNets.cs | 300 ++++++++++++++++++ .../GameCult.Geometry.csproj | 1 + .../SURFACE-NETS-PROVENANCE.md | 36 +++ 5 files changed, 490 insertions(+) create mode 100644 src/GameCult.Geometry/CultGeometryQuadMesh.cs create mode 100644 src/GameCult.Geometry/CultGeometryQuadNormals.cs create mode 100644 src/GameCult.Geometry/CultGeometrySurfaceNets.cs create mode 100644 src/GameCult.Geometry/SURFACE-NETS-PROVENANCE.md diff --git a/src/GameCult.Geometry/CultGeometryQuadMesh.cs b/src/GameCult.Geometry/CultGeometryQuadMesh.cs new file mode 100644 index 00000000..ab40ba4a --- /dev/null +++ b/src/GameCult.Geometry/CultGeometryQuadMesh.cs @@ -0,0 +1,88 @@ +using System; + +namespace GameCult.Geometry +{ + /// + /// The grid edge from sample (X, Y, Z) to (X, Y, Z) + e_Axis, where Axis 0/1/2 is X/Y/Z. + /// + public readonly struct CultGeometryGridEdge : IEquatable, IComparable + { + /// Creates a grid edge. + public CultGeometryGridEdge(int x, int y, int z, byte axis) + { + X = x; + Y = y; + Z = z; + Axis = axis; + } + + /// Gets the starting sample's x coordinate. + public int X { get; } + + /// Gets the starting sample's y coordinate. + public int Y { get; } + + /// Gets the starting sample's z coordinate. + public int Z { get; } + + /// Gets the axis the edge runs along: 0 for x, 1 for y, 2 for z. + public byte Axis { get; } + + /// Returns whether this edge is equal to another edge. + public bool Equals(CultGeometryGridEdge other) => + X == other.X && Y == other.Y && Z == other.Z && Axis == other.Axis; + + /// + public override bool Equals(object? obj) => obj is CultGeometryGridEdge other && Equals(other); + + /// + public override int GetHashCode() => HashCode.Combine(X, Y, Z, Axis); + + /// Compares edges by (X, Y, Z, Axis) in that order. + public int CompareTo(CultGeometryGridEdge other) + { + var byX = X.CompareTo(other.X); + if (byX != 0) return byX; + + var byY = Y.CompareTo(other.Y); + if (byY != 0) return byY; + + var byZ = Z.CompareTo(other.Z); + if (byZ != 0) return byZ; + + return Axis.CompareTo(other.Axis); + } + + /// + public override string ToString() => $"({X}, {Y}, {Z})+e{Axis}"; + + public static bool operator ==(CultGeometryGridEdge left, CultGeometryGridEdge right) => left.Equals(right); + public static bool operator !=(CultGeometryGridEdge left, CultGeometryGridEdge right) => !left.Equals(right); + public static bool operator <(CultGeometryGridEdge left, CultGeometryGridEdge right) => left.CompareTo(right) < 0; + public static bool operator >(CultGeometryGridEdge left, CultGeometryGridEdge right) => left.CompareTo(right) > 0; + public static bool operator <=(CultGeometryGridEdge left, CultGeometryGridEdge right) => left.CompareTo(right) <= 0; + public static bool operator >=(CultGeometryGridEdge left, CultGeometryGridEdge right) => left.CompareTo(right) >= 0; + } + + /// + /// An engine-neutral quad mesh with welded vertices, produced by surface nets extraction. + /// Not a CultCache document: this is an in-process algorithm result, not a persisted artifact. + /// + public sealed class CultGeometryQuadMesh + { + /// Gets or sets the welded vertex positions, xyz per vertex. + public float[] Positions { get; set; } = Array.Empty(); + + /// Gets or sets the quad vertex indices, 4 per quad. + public uint[] Quads { get; set; } = Array.Empty(); + + /// Gets or sets the grid edge each quad was emitted for, one per quad, strictly ascending. + public CultGeometryGridEdge[] QuadEdges { get; set; } = Array.Empty(); + + /// Gets the number of quads. + public int QuadCount => Quads.Length / 4; + + /// Gets the number of welded vertices. + public int VertexCount => Positions.Length / 3; + } +} diff --git a/src/GameCult.Geometry/CultGeometryQuadNormals.cs b/src/GameCult.Geometry/CultGeometryQuadNormals.cs new file mode 100644 index 00000000..a8121ca7 --- /dev/null +++ b/src/GameCult.Geometry/CultGeometryQuadNormals.cs @@ -0,0 +1,65 @@ +using System; +using CultMath; + +namespace GameCult.Geometry +{ + /// + /// Derives per-vertex normals for a welded quad mesh by face-weighted averaging. + /// + public static class CultGeometryQuadNormals + { + /// + /// Computes a per-vertex normal for as the sum, over every quad + /// touching that vertex, of the quad's diagonal cross product (which points along the + /// quad's normal and whose length is twice a non-planar quad's area), then normalized. + /// A vertex touched only by zero-area quads gets a zero normal. + /// + public static float[] FaceWeighted(CultGeometryQuadMesh mesh) + { + if (mesh == null) throw new ArgumentNullException(nameof(mesh)); + + var sums = new float3[mesh.VertexCount]; + for (var quad = 0; quad < mesh.QuadCount; quad++) + { + var baseIndex = quad * 4; + var i0 = mesh.Quads[baseIndex]; + var i1 = mesh.Quads[baseIndex + 1]; + var i2 = mesh.Quads[baseIndex + 2]; + var i3 = mesh.Quads[baseIndex + 3]; + + var p0 = Position(mesh.Positions, i0); + var p1 = Position(mesh.Positions, i1); + var p2 = Position(mesh.Positions, i2); + var p3 = Position(mesh.Positions, i3); + + var faceVector = math.cross(p2 - p0, p3 - p1); + + sums[i0] += faceVector; + sums[i1] += faceVector; + sums[i2] += faceVector; + sums[i3] += faceVector; + } + + var normals = new float[mesh.VertexCount * 3]; + for (var vertex = 0; vertex < sums.Length; vertex++) + { + var sum = sums[vertex]; + + // A vertex whose contributions cancel or were all zero-area gets a zero normal; + // normalize follows HLSL and would return NaN for it (matches Extract's degenerate rule). + var normal = math.dot(sum, sum) > 0f ? math.normalize(sum) : float3.zero; + normals[vertex * 3] = normal.x; + normals[vertex * 3 + 1] = normal.y; + normals[vertex * 3 + 2] = normal.z; + } + + return normals; + } + + private static float3 Position(float[] positions, uint index) + { + var offset = (int)index * 3; + return new float3(positions[offset], positions[offset + 1], positions[offset + 2]); + } + } +} diff --git a/src/GameCult.Geometry/CultGeometrySurfaceNets.cs b/src/GameCult.Geometry/CultGeometrySurfaceNets.cs new file mode 100644 index 00000000..c7df880f --- /dev/null +++ b/src/GameCult.Geometry/CultGeometrySurfaceNets.cs @@ -0,0 +1,300 @@ +using System; +using System.Collections.Generic; +using CultMath; + +namespace GameCult.Geometry +{ + /// + /// Extracts engine-neutral welded quad meshes from sampled scalar fields with surface nets. + /// + public static class CultGeometrySurfaceNets + { + // Start-corner offsets for each axis's four parallel edges of a unit cell. The edge itself + // runs from the listed corner to that corner plus one step along Axis. + private static readonly (int Dx, int Dy, int Dz)[] AxisXEdgeStarts = + { + (0, 0, 0), (0, 1, 0), (0, 0, 1), (0, 1, 1), + }; + + private static readonly (int Dx, int Dy, int Dz)[] AxisYEdgeStarts = + { + (0, 0, 0), (1, 0, 0), (0, 0, 1), (1, 0, 1), + }; + + private static readonly (int Dx, int Dy, int Dz)[] AxisZEdgeStarts = + { + (0, 0, 0), (1, 0, 0), (0, 1, 0), (1, 1, 0), + }; + + /// + /// Extracts the surface with surface nets. + /// Values less than or equal to the isovalue are treated as inside. + /// + public static CultGeometryQuadMesh Extract( + float[,,] samples, + float isoValue = 0f, + CultVec3 origin = default, + float cellSize = 1f) + { + if (samples == null) throw new ArgumentNullException(nameof(samples)); + if (samples.GetLength(0) < 2 || samples.GetLength(1) < 2 || samples.GetLength(2) < 2) + { + throw new ArgumentException("A surface nets field requires at least two samples on every axis.", nameof(samples)); + } + + if (!(cellSize > 0f) || float.IsInfinity(cellSize)) + { + throw new ArgumentOutOfRangeException(nameof(cellSize), "Cell size must be finite and positive."); + } + + var sizeX = samples.GetLength(0); + var sizeY = samples.GetLength(1); + var sizeZ = samples.GetLength(2); + var cellsX = sizeX - 1; + var cellsY = sizeY - 1; + var cellsZ = sizeZ - 1; + var originValue = new float3(origin.X, origin.Y, origin.Z); + + var vertexIndex = new Dictionary<(int, int, int), uint>(); + var positions = new List(); + + for (var cx = 0; cx < cellsX; cx++) + for (var cy = 0; cy < cellsY; cy++) + for (var cz = 0; cz < cellsZ; cz++) + { + var vertex = CellVertex(samples, cx, cy, cz, isoValue, originValue, cellSize); + if (vertex is not { } position) continue; + + vertexIndex[(cx, cy, cz)] = (uint)(positions.Count / 3); + Append(position, positions); + } + + var quads = new List(); + var quadEdges = new List(); + + for (var x = 0; x < sizeX; x++) + for (var y = 0; y < sizeY; y++) + for (var z = 0; z < sizeZ; z++) + { + for (byte axis = 0; axis < 3; axis++) + { + if (!TryGetOtherEndpoint(sizeX, sizeY, sizeZ, x, y, z, axis, out var ox, out var oy, out var oz)) + { + continue; + } + + var inside = samples[x, y, z] <= isoValue; + var otherInside = samples[ox, oy, oz] <= isoValue; + if (inside == otherInside) continue; + + if (!TryGetQuadCells(cellsX, cellsY, cellsZ, x, y, z, axis, out var c0, out var c1, out var c2, out var c3)) + { + continue; + } + + if (!vertexIndex.TryGetValue(c0, out var v0) || + !vertexIndex.TryGetValue(c1, out var v1) || + !vertexIndex.TryGetValue(c2, out var v2) || + !vertexIndex.TryGetValue(c3, out var v3)) + { + // The four cells around a crossing edge each contain that same crossing edge + // among their own twelve, so each of them is mixed and has a vertex. This is an + // invariant guard, not a reachable branch. + throw new InvalidOperationException("A surface-nets crossing edge bordered a cell with no vertex."); + } + + var p0 = Position(positions, v0); + var p1 = Position(positions, v1); + var p2 = Position(positions, v2); + var p3 = Position(positions, v3); + var diagonalCross = math.cross(p2 - p0, p3 - p1); + + var insideCoord = inside ? new float3(x, y, z) : new float3(ox, oy, oz); + var outsideCoord = inside ? new float3(ox, oy, oz) : new float3(x, y, z); + var outward = outsideCoord - insideCoord; + + if (math.dot(diagonalCross, outward) < 0f) + { + (v1, v3) = (v3, v1); + } + + quads.Add(v0); + quads.Add(v1); + quads.Add(v2); + quads.Add(v3); + quadEdges.Add(new CultGeometryGridEdge(x, y, z, axis)); + } + } + + return new CultGeometryQuadMesh + { + Positions = positions.ToArray(), + Quads = quads.ToArray(), + QuadEdges = quadEdges.ToArray(), + }; + } + + private static float3? CellVertex( + float[,,] samples, + int cx, + int cy, + int cz, + float isoValue, + float3 origin, + float cellSize) + { + var sum = float3.zero; + var count = 0; + + AccumulateAxisEdges(samples, cx, cy, cz, 0, AxisXEdgeStarts, isoValue, origin, cellSize, ref sum, ref count); + AccumulateAxisEdges(samples, cx, cy, cz, 1, AxisYEdgeStarts, isoValue, origin, cellSize, ref sum, ref count); + AccumulateAxisEdges(samples, cx, cy, cz, 2, AxisZEdgeStarts, isoValue, origin, cellSize, ref sum, ref count); + + return count == 0 ? null : sum / count; + } + + private static void AccumulateAxisEdges( + float[,,] samples, + int cx, + int cy, + int cz, + int axis, + (int Dx, int Dy, int Dz)[] starts, + float isoValue, + float3 origin, + float cellSize, + ref float3 sum, + ref int count) + { + foreach (var start in starts) + { + var startX = cx + start.Dx; + var startY = cy + start.Dy; + var startZ = cz + start.Dz; + var endX = startX + (axis == 0 ? 1 : 0); + var endY = startY + (axis == 1 ? 1 : 0); + var endZ = startZ + (axis == 2 ? 1 : 0); + + var startValue = samples[startX, startY, startZ]; + var endValue = samples[endX, endY, endZ]; + if ((startValue <= isoValue) == (endValue <= isoValue)) continue; + + var startPosition = origin + new float3(startX, startY, startZ) * cellSize; + var endPosition = origin + new float3(endX, endY, endZ) * cellSize; + sum += Interpolate(startPosition, startValue, endPosition, endValue, isoValue); + count++; + } + } + + private static float3 Interpolate(float3 first, float firstValue, float3 second, float secondValue, float isoValue) + { + var delta = secondValue - firstValue; + var amount = math.abs(delta) <= 1e-20f ? 0.5f : (isoValue - firstValue) / delta; + return math.lerp(first, second, amount); + } + + private static bool TryGetOtherEndpoint( + int sizeX, + int sizeY, + int sizeZ, + int x, + int y, + int z, + byte axis, + out int ox, + out int oy, + out int oz) + { + ox = x; + oy = y; + oz = z; + + switch (axis) + { + case 0: + if (x + 1 >= sizeX) return false; + ox = x + 1; + return true; + case 1: + if (y + 1 >= sizeY) return false; + oy = y + 1; + return true; + default: + if (z + 1 >= sizeZ) return false; + oz = z + 1; + return true; + } + } + + // The four cells that share the grid edge from (x, y, z) to (x, y, z) + e_axis, listed as a + // perimeter loop around the edge (not crossed diagonally) so the caller can build a quad + // directly from them. + private static bool TryGetQuadCells( + int cellsX, + int cellsY, + int cellsZ, + int x, + int y, + int z, + byte axis, + out (int, int, int) c0, + out (int, int, int) c1, + out (int, int, int) c2, + out (int, int, int) c3) + { + c0 = c1 = c2 = c3 = default; + + switch (axis) + { + case 0: + if (!InRange(y - 1, cellsY) || !InRange(y, cellsY) || !InRange(z - 1, cellsZ) || !InRange(z, cellsZ)) + { + return false; + } + + c0 = (x, y - 1, z - 1); + c1 = (x, y, z - 1); + c2 = (x, y, z); + c3 = (x, y - 1, z); + return InRange(x, cellsX); + case 1: + if (!InRange(x - 1, cellsX) || !InRange(x, cellsX) || !InRange(z - 1, cellsZ) || !InRange(z, cellsZ)) + { + return false; + } + + c0 = (x - 1, y, z - 1); + c1 = (x, y, z - 1); + c2 = (x, y, z); + c3 = (x - 1, y, z); + return InRange(y, cellsY); + default: + if (!InRange(x - 1, cellsX) || !InRange(x, cellsX) || !InRange(y - 1, cellsY) || !InRange(y, cellsY)) + { + return false; + } + + c0 = (x - 1, y - 1, z); + c1 = (x, y - 1, z); + c2 = (x, y, z); + c3 = (x - 1, y, z); + return InRange(z, cellsZ); + } + } + + private static bool InRange(int value, int count) => value >= 0 && value < count; + + private static float3 Position(List positions, uint index) + { + var offset = (int)index * 3; + return new float3(positions[offset], positions[offset + 1], positions[offset + 2]); + } + + private static void Append(float3 value, List destination) + { + destination.Add(value.x); + destination.Add(value.y); + destination.Add(value.z); + } + } +} diff --git a/src/GameCult.Geometry/GameCult.Geometry.csproj b/src/GameCult.Geometry/GameCult.Geometry.csproj index a3ad9457..c10c2ecf 100644 --- a/src/GameCult.Geometry/GameCult.Geometry.csproj +++ b/src/GameCult.Geometry/GameCult.Geometry.csproj @@ -17,5 +17,6 @@ + diff --git a/src/GameCult.Geometry/SURFACE-NETS-PROVENANCE.md b/src/GameCult.Geometry/SURFACE-NETS-PROVENANCE.md new file mode 100644 index 00000000..989f9ae9 --- /dev/null +++ b/src/GameCult.Geometry/SURFACE-NETS-PROVENANCE.md @@ -0,0 +1,36 @@ +# Surface nets extraction provenance + +`CultGeometrySurfaceNets` is a clean implementation of surface nets over a +regular scalar grid. It is implemented fresh from the published algorithm +description; no code or lookup tables were copied from any source. `CultGeometryQuadNormals` +is likewise a fresh implementation of face-weighted vertex normal averaging. + +Cited literature: + +- S. F. F. Gibson, "Constrained Elastic Surface Nets: Generating Smooth + Surfaces from Binary Segmented Data," MICCAI 1998 (MERL Technical Report + TR99-24). +- M. Lysenko, "Smooth Voxel Terrain (Part 2)," 0fps.net, 2012. +- N. Max, "Weights for Computing Vertex Normals from Facet Normals," Journal + of Graphics Tools 4(2), 1999. + +Deliberate behavior: + +- a scalar sample less than or equal to the isovalue is inside, matching + `CultGeometryIsoSurface`'s sign convention; +- a cell's vertex is the mean of the linearly interpolated crossings of its + own twelve edges; +- Gibson's elastic relaxation of the vertex toward its neighbors is omitted; + the vertex stays at the unrelaxed dual-cell mean; +- one quad is emitted per crossing interior grid edge, from the (up to) four + cells that share it; an edge whose four surrounding cells do not all exist + emits no quad, leaving the mesh open at the field boundary; +- quad winding is oriented outward: the diagonal cross product of each quad + must have a positive dot product with the direction from the edge's inside + endpoint to its outside endpoint; +- vertices are welded: one vertex per active cell, shared by every quad + touching that cell; +- face-weighted normals sum each quad's diagonal cross product (unnormalized, + so a bigger quad outweighs a sliver) per touched vertex, then normalize; a + vertex touched only by zero-area quads gets a zero normal rather than NaN, + matching `CultGeometryIsoSurface`'s degenerate-normal rule.