SpatialMapping.Core 0.3.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.

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.

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.

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.