SpatialMapping.Core 0.2.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);

// 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 others (continuous-ray-vs-point/segment hits are measure-zero in 3D — caller should fall back to AABB-based FindFirstRayHit).

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.