SpatialMapping.Core 0.5.0

SpatialMapping.Core

Spatial indices and geometry primitives for 3D (and 2D) work in .NET. Six index structures, a unified query set (nearest / k-NN / radius / segment / ray / box / all-pairs), a full pure-geometry layer in both float and double precision (planes, triangles, spheres, polylines, rigid transforms, best-fit), and an empirical benchmark harness that measures every combination on your data so you can pick by reading numbers, not guessing.

Multi-targets netstandard2.0 + net8.0. Works in modern .NET, .NET Framework 4.7.2+, Mono, Unity.

Why this library, not [X]?

The C# spatial-geometry landscape is fragmented. SpatialMapping.Core sits in the empty middle: not GIS-oriented like NetTopologySuite, not mesh-oriented like geometry3sharp, not primitives-only like Math.NET Spatial, not single-structure like the various k-d / R-tree libraries. The point is to answer "I have N items in 3D — what do I do with them?" in one package, with one API, and let an empirical benchmark pick the right structure.

Pick a different library if you need:

Pick this library if you need:

• Multi-structure spatial indices under one API. Switch between BVH, k-d tree, octree, AABB grid, uniform grid, and brute-force without rewriting query code.
• Empirical "which structure is fastest on my data" measurement. SpatialBenchmark.RunForSegments(yourSegments) runs every applicable index × every query type and returns a markdown table you can paste into a PR.
• Float + double parity across the pure-geometry layer. Most alternatives are float-only or double-only; this package ships both with explicit conversion APIs.
• netstandard2.0 reach. Modern .NET, but also .NET Framework 4.7.2 (Revit / AutoCAD addins), Mono, Unity. The netstandard2.0 binary is runtime-tested as part of the test suite, not just compile-checked.
• Tight scope. No GIS DSL, no engine dependency, no mesh-CSG ambition. Fewer concepts to learn; one nupkg with no transitive dependencies on netstandard2.0 (and only System.Numerics.Vectors + System.Memory on the netstandard2.0 target).

If you'd reach for both NetTopologySuite and geometry3sharp on the same project, you're not the target user. If you'd reach for neither — just want spatial indices that work — this is for you.

Quick start

Points

using SpatialMapping.Core;
using System.Numerics;

var positions = new[] { new Vector3(0,0,0), new Vector3(1,0,0), new Vector3(0,5,0) };
var bvh = new BvhIndex3D<Vector3, Vector3Distance>();
bvh.Build(positions);

var nearest = bvh.FindNearestToPoint(new Vector3(0.4f, 0, 0));   // .Index = 0

var topK = new List<NearestResult<Vector3>>();
bvh.FindKNearestToPoint(new Vector3(0,0,0), k: 2, topK);

var withinR = new List<NearestResult<Vector3>>();
bvh.FindWithinRadiusOfPoint(Vector3.Zero, radius: 3f, withinR);

Segments (or any AABB-shaped item)

var segments = new[]
{
    new Segment3D(new Vector3(0,0,0), new Vector3(1000, 0, 0)), // long
    new Segment3D(new Vector3(0,5,0), new Vector3(0, 5, 1)),    // short
};
var bvh = new BvhIndex3D<Segment3D, Segment3DDistance>();
bvh.Build(segments);

// Endpoints are 300 and 700 units from this query, but distance to the segment is 0.
var hit = bvh.FindNearestToPoint(new Vector3(300, 0, 0));        // .Index = 0, d² = 0

var clash = bvh.FindNearestToSegment(new Segment3D(new Vector3(0, -1, 0), new Vector3(0, 1, 0)));

var inBox = new List<NearestResult<Segment3D>>();
bvh.QueryBox(new Bounds3D(new Vector3(-1,-1,-1), new Vector3(2,2,2)), inBox);

var firstHit = bvh.FindFirstRayHit(new Vector3(-5, 0, 0), new Vector3(1, 0, 0));

All-pairs queries (clash detection, collision broad-phase, self-intersection)

// One call instead of N radius queries with HashSet-dedup.
var pairs = new List<PairResult<Segment3D>>();
bvh.FindAllPairsWithin(radius: 0.05f, pairs);
// Each unordered pair appears exactly once with pairs[i].IndexA < pairs[i].IndexB.

// Cross-index spatial join: which beams clash with which columns?
var beamColumnPairs = new List<PairResult<Segment3D>>();
beamBvh.FindAllPairsWithin(columnBvh, radius: 0.05f, beamColumnPairs);

Tight ray-vs-item (genuine geometric hits, not AABB hits)

// Triangles support Möller-Trumbore ray intersection out of the box.
var tris = new BvhIndex3D<Triangle3D, Triangle3DDistance>();
tris.Build(meshTriangles);
var hit = tris.FindFirstRayHitTight(rayOrigin, rayDirection);
// hit.Distance is the parametric t at the real triangle intersection.

// Spheres solve the ray-sphere quadratic.
var balls = new BvhIndex3D<Sphere3D, Sphere3DDistance>();
balls.Build(spheres);
var nearestSphere = balls.FindFirstRayHitTight(rayOrigin, rayDirection);

Empirical "which structure?" — SpatialBenchmark

[Fact]
public void PickIndexForMyData()
{
    var report = SpatialBenchmark.RunForSegments(myActualData);
    _out.WriteLine(report.ToMarkdown());
    Assert.Empty(report.CorrectnessIssues);
}

The harness times every applicable structure on every supported query. Naive O(N) fallbacks fill in queries a structure doesn't natively support (so every cell has a number — the timing tells you the cost of forcing the wrong tool). A per-query-type timeout lets you cap infeasible combinations.

Geometry primitives

A full narrow-phase pure-geometry layer ships alongside the indices. Both float (default) and double-precision (*d suffix) types are provided for callers who need machine precision in narrow-phase math (e.g., converting between coordinate systems with large translations, geometric cleanup).

3D primitives

// Plane: factories cover the three common ways to build one.
var p1 = Plane3D.FromPointAndNormal(new Vector3(0, 0, 10), Vector3.UnitZ);
var p2 = Plane3D.FromThreePoints(a, b, c);
var p3 = Plane3D.FromPointAndTwoVectors(origin, ux, uy);

float signed = p1.SignedDistanceTo(point);
Vector3 projection = p1.ClosestPointTo(point);
Vector3 mirrored = p1.Reflect(point);

// Triangle: closest-point with barycentric (Ericson §5.1.5).
var t = new Triangle3D(v0, v1, v2);
var closest = t.ClosestPointTo(p, out float u, out float v, out float w);
// u + v + w == 1, point = u*v0 + v*v1 + w*v2.

// Sphere, oriented bounding box, polyline (open or closed).
var s = new Sphere3D(center, radius);
var obb = OrientedBounds3D.FromPoints(pointSet);            // PCA-fit via Jacobi
var pl = new Polyline3D(vertices, isClosed: true);
// Closed polylines expose a wrap-around segment: `pl.SegmentCount == vertices.Count`
// (one more than the open case), and `pl.GetSegment(pl.SegmentCount - 1)` connects
// the last vertex back to the first.

// Intersection routines: line-plane, plane-plane, three-plane, segment-plane, segment-triangle (Möller-Trumbore).
if (Intersect.SegmentWithTriangle(seg, tri, out float t_, out Vector3 point, out float bu, out float bv, out float bw)) { /* ... */ }

Rigid transforms

// Build a coordinate-frame conversion once, apply it to every primitive.
var revitToSap = Transform3D.FromBasis(origin, axisX, axisY, axisZ)
               * Transform3D.UniformScale(0.3048f);  // feet → meters

var pSap     = revitToSap.Apply(pRevit);              // point
var dSap     = revitToSap.ApplyDirection(dRevit);     // direction (no translation)
var segSap   = revitToSap.Apply(segRevit);            // Segment3D
var planeSap = revitToSap.Apply(planeRevit);          // Plane3D
var triSap   = revitToSap.Apply(triRevit);            // Triangle3D
var aabbSap  = revitToSap.Apply(aabbRevit);           // Bounds3D (eight-corner trick)
var inverse  = revitToSap.Inverse();

// Predefined factories.
Transform3D.YUpToZUp();
Transform3D.LookAt(eye, target, up);
Transform3D.FromAxisAngle(Vector3.UnitZ, MathF.PI / 4);

Quaternion is System.Numerics.Quaternion (float); Quaterniond is the double-precision sibling with slerp, from-axis-angle, from-two-vectors, rotation-matrix conversion, and the standard Rotate(v) = q · v · q⁻¹ operation. QuaternionExtensions.FromTwoVectors / FromRotationMatrix fills in the gaps in System.Numerics.

Best-fit & PCA

// Line through points: largest-eigenvalue eigenvector of the covariance.
var line = BestFit.Line(points);
// line.Origin, line.UnitDirection, line.MaxResidual, line.RmsResidual.

// Plane through points: smallest-eigenvalue eigenvector.
var plane = BestFit.Plane(points);

// Sphere (algebraic 4×4 normal equations) and circle (plane fit + 2D algebraic).
var sphereFit = BestFit.Sphere(points);
var circleFit = BestFit.Circle(coplanarPoints);

// Full eigendecomposition (centroid + major/mid/minor axes + eigenvalues).
var axes = BestFit.PrincipalAxes(points);

Snap / weld / quantize (machine-precision cleanup)

// Round each coordinate to the nearest multiple of a grid.
Vector3 snapped = Snap.Quantize(point, gridSize: 0.001f);   // 1mm grid

// Cluster near-coincident points within tolerance, collapse to centroid.
Vector3[] welded = Snap.WeldPoints(pointSet, tolerance: 0.005f);

// Same operation on segment endpoints — produces FEA-ready joint connectivity.
Segment3D[] joinedUp = Snap.WeldEndpoints(segments, tolerance: 0.01f);

Clustering is transitive: if A is within ε of B and B is within ε of C, all three merge — even when A and C are slightly farther apart than ε. Implementation uses a BVH for broad-phase candidate finding plus union-find for cluster construction.

2D layer

Vector2d, Bounds2D / Bounds2Dd, Segment2D / Segment2Dd, BezierCurve2D / BezierCurve2Dd (with adaptive chord-deviation tessellation), Polyline2D / Polyline2Dd, ConvexHull2D (Andrew's monotone chain), Delaunay2D.Triangulate (Bowyer-Watson), and BentleyOttmann2D.FindIntersections (sweep-line all-pairs segment intersection in O((N + K) log N)).

A heuristic PathClassifier takes already-extracted 2D path geometry (line segments + cubic Béziers + closed flag — format-agnostic, supplied by your own PDF/SVG/DXF reader) and tags each subpath as Triangle, Rectangle, RevisionCloud, Polygon, Polyline, or Empty with tunable thresholds.

Mesh queries

// Closed-mesh inside/outside test via ray-parity through a triangle BVH.
var triBvh = new BvhIndex3D<Triangle3D, Triangle3DDistance>();
triBvh.Build(closedMeshTriangles);
bool inside = Mesh.IsInside(triBvh, point);

Pick a structure

Empirically: at small N (~few hundred), AabbGridIndex3D can beat BvhIndex3D on total time because its near-constant-time build (one pass vs O(N log N) sort) outweighs query-time differences. The Query × Structure matrix below tells you about per-query capabilities; the SpatialBenchmark harness tells you what's actually fastest on your data.

If unsure, run SpatialBenchmark on your data and read the table.

Query × structure matrix

Native = the structure's own implementation. Naive = O(N) fallback (still correct, just slow). "—" = not supported by this structure (use a different one).

All six expose batch + parallel variants of FindNearest* for read-only-after-build query workloads.

Adapting custom item types

Implement IItemDistance3D<T> as a readonly struct for the JIT-monomorphized fast path:

public readonly struct SphereDistance : IItemDistance3D<Sphere>
{
    public Bounds3D GetAabb(in Sphere s)
        => new(s.Center - new Vector3(s.Radius), s.Center + new Vector3(s.Radius));

    public float DistanceSquaredToPoint(in Sphere s, Vector3 p)
    {
        float d = Vector3.Distance(s.Center, p) - s.Radius;
        return d <= 0 ? 0 : d * d;
    }

    public float DistanceSquaredToSegment(in Sphere s, in Segment3D q)
    {
        float d = MathF.Sqrt(q.DistanceSquaredTo(s.Center)) - s.Radius;
        return d <= 0 ? 0 : d * d;
    }

    public float DistanceSquaredToItem(in Sphere a, in Sphere b)
    {
        // Used by all-pairs / dual-BVH queries.
        float d = Vector3.Distance(a.Center, b.Center) - a.Radius - b.Radius;
        return d <= 0f ? 0f : d * d;
    }

    public bool TryRayHit(in Sphere s, Vector3 o, Vector3 dir, float maxDistance, out float t)
    {
        // Used by FindFirstRayHitTight / FindAllRayHitsTight. Return false for item
        // types where "tight ray hit" isn't meaningful (e.g., 3D points).
        // ... ray-sphere quadratic ...
    }
}

var bvh = new BvhIndex3D<Sphere, SphereDistance>();
bvh.Build(spheres);

Or use the delegate-flavored convenience overload (BvhIndex3D<Sphere>(getAabb, distToPoint, distToSegment) plus optional distToItem and tryRayHit 4- and 5-arg constructors) for ad-hoc work — one extra delegate dispatch per leaf hit.

Built-in adapters: Vector3Distance, Segment3DDistance, Sphere3DDistance, Triangle3DDistance, Polyline3DDistance. Pair queries are supported on all of them; tight ray hit is supported on Triangle3DDistance (Möller-Trumbore) and Sphere3DDistance (quadratic), and is a no-op on the other three — Vector3Distance, Segment3DDistance, Polyline3DDistance. Points, segments, and polylines are all 0- or 1-dimensional in 3D, so a continuous ray hits them with probability zero — there's nothing meaningful to report. Use AABB-based FindFirstRayHit instead, or build a thick-ray query on top of the adapter's DistanceSquaredToSegment.

The "Tight ray" column in the matrix above means the BVH supports the dispatch — whether it returns hits is a property of the adapter, not the structure.

Factory conventions

The package uses five construction patterns. The naming reflects whether the factory derives (does math), carries a precondition (named contract on inputs), matches the field shape of the underlying record, is a canonical constant, or is a canonical zero-arg method.

Quick decision tree:

• "I have raw inputs and want the library to derive the math" → FromXxx derivation.
• "I have field-shaped inputs but the library has a named precondition" → FromXxx precondition.
• "I have field-shaped inputs with no special precondition" → Create.
• "I want a canonical instance" → static property if it's a constant, bare-verb method if it allocates.

Plane3D and Transform3D carry both derivation-FromXxx and field-shape factories because both add value: derivation handles raw input, field-shape lets callers who've done the math themselves skip the redundant normalization step. The "Naming convention for constructor named arguments" subsection below explains why the field-shape factories exist at all (records force PascalCase parameter names; the factories give camelCase).

Naming convention for constructor named arguments

Most geometry primitives in this package are readonly record structs with a positional constructor. C# records inherit the property casing (PascalCase) for parameter names, not the usual camelCase. So when you reach for a named argument, type it PascalCase:

// works
var s = new Sphere3D(Center: new Vector3(0, 0, 0), Radius: 10f);
var t = new Triangle3D(V0: a, V1: b, V2: c);
var p = new Plane3D(Normal: n, D: 5f);

// CS1739: the parameter is named "Radius", not "radius"
var s2 = new Sphere3D(center: new Vector3(0, 0, 0), radius: 10f);

Affects every positional record in this package — Sphere3D, Segment3D, Triangle3D, Plane3D, Bounds3D, Circle3D, OrientedBounds3D, Sphere3Dd, … and the 2D / path-classifier pieces LinePiece2D, BezierPiece2D. IntelliSense will show the correct names in the parameter hint; the lower-case form just won't compile.

A note on in parameters

Methods that take a struct of more than two Vector3s (e.g. Segment3D, Plane3D, Triangle3D, Bounds3D) declare those parameters as in to avoid copying on call. Writing in at the call site is optional — the C# compiler will pass by readonly reference either way:

var seg = new Segment3D(a, b);
sphere.DistanceSquaredTo(seg);     // fine — `in` is implicit
sphere.DistanceSquaredTo(in seg);  // also fine, just explicit

So the in you see in IntelliSense doesn't force any extra syntax at the call site.

Helper classes amenable to using static

Most calculation helpers in this package are stateless static classes — pull them in with using static at file scope when the call sites are noisy:

using static SpatialMapping.Core.MathHelpers;
using static SpatialMapping.Core.Snap;
using static SpatialMapping.Core.SegmentMath;

var z = Clamp(value, 0, 100);          // MathHelpers.Clamp
var snapped = Quantize(p, 0.001f);     // Snap.Quantize
var d2 = SquaredDistance(segA, segB);  // SegmentMath.SquaredDistance

The full list:

The pure-geometry primitives (Plane3D, Triangle3D, Sphere3D, Transform3D, …) are types rather than helper classes — call their methods on the value or use the static FromXxx / Create factories.

Thread safety

After a successful Build (or Refit), a spatial index is read-only and safe to query concurrently from any number of threads. The traversal uses only stack-local state and never mutates the tree. So this pattern is fine:

var bvh = new BvhIndex3D<Segment3D, Segment3DDistance>();
bvh.Build(segments);
// no further mutation of bvh anywhere

Parallel.For(0, queries.Count, i =>
{
    var hit = bvh.FindNearestToPoint(queries[i]);          // safe
    // …record hit somewhere thread-local…
});

What is not safe across threads:

• Concurrent Build and queries on the same instance. Build mutates the internal arrays. Don't query an instance while another thread is still building it. (Building multiple separate instances on different threads is fine.)
• Sharing a List<T> results output across threads. Methods like FindWithinRadiusOfPoint(query, radius, results) clear and fill results — if two threads pass the same list, they'll race on it. Give each thread its own list, or use the array-backed …Parallel variants (FindNearestToPointBatchParallel, FindNearestToSegmentBatchParallel) where each thread writes into its own slot.
• Concurrent Refit and queries. Same as Build — Refit mutates AABB caches.

The pure-geometry layer (Plane3D, Triangle3D, Transform3D, Quaterniond, etc.) is value types and free functions — no shared mutable state, threadsafe by default.

Out of scope (today)

• No serialization. In-process use only; rebuild from source data on load.
• No frustum query. Same traversal as ray; trivial to add when needed.
• No parallel build. Rebuild is single-threaded; query is parallel.
• Indices are float-only (System.Numerics.Vector3). The pure-geometry layer offers double-precision mirrors (Vector3d, Segment3Dd, Plane3Dd, etc.) for narrow-phase math that needs machine precision.
• No 3D convex hull, best-fit cylinder, or marching cubes yet — flagged in the design notes as the next incremental additions.

Versioning

0.x — early. API may change between minor versions. The 0.1 → 0.2 transition added the post-build all-pairs primitive and the entire pure-geometry layer, with two non-backwards-compatible extensions to IItemDistance3D<T> (DistanceSquaredToItem and TryRayHit); custom-adapter implementers need to add the new methods. Built-in adapters were updated.