From 2214016e6f5ceed8228692c4e8ebfc738e6cf022 Mon Sep 17 00:00:00 2001 From: Jonas Kamsker <11245306+JKamsker@users.noreply.github.com> Date: Sun, 28 Sep 2025 03:31:59 +0200 Subject: [PATCH 1/8] Add spatial data structures and functionality - Implemented GeoBoundingBox for spatial bounding box calculations. - Added GeoJson for serialization and deserialization of GeoShape objects. - Introduced GeoMath for geographical calculations including distance metrics. - Created GeoShape as an abstract base class with derived types: GeoPoint, GeoLineString, and GeoPolygon. - Added GeoValidation for coordinate validation. - Developed Geometry utility class for geometric operations like point containment and intersection checks. - Implemented LongitudeRange for handling longitude ranges with wrapping. - Created Spatial class for spatial indexing and querying capabilities. - Added SpatialIndexing for Morton code computation for spatial indexing. - Introduced SpatialOptions for configuration of spatial operations. --- LiteDB.Tests/Spatial/SpatialTests.cs | 197 ++++++++++++++ LiteDB/Spatial/GeoBoundingBox.cs | 77 ++++++ LiteDB/Spatial/GeoJson.cs | 194 ++++++++++++++ LiteDB/Spatial/GeoMath.cs | 126 +++++++++ LiteDB/Spatial/GeoShape.cs | 163 ++++++++++++ LiteDB/Spatial/GeoValidation.cs | 25 ++ LiteDB/Spatial/Geometry.cs | 293 +++++++++++++++++++++ LiteDB/Spatial/LongitudeRange.cs | 67 +++++ LiteDB/Spatial/Spatial.cs | 378 +++++++++++++++++++++++++++ LiteDB/Spatial/SpatialIndexing.cs | 52 ++++ LiteDB/Spatial/SpatialOptions.cs | 25 ++ 11 files changed, 1597 insertions(+) create mode 100644 LiteDB.Tests/Spatial/SpatialTests.cs create mode 100644 LiteDB/Spatial/GeoBoundingBox.cs create mode 100644 LiteDB/Spatial/GeoJson.cs create mode 100644 LiteDB/Spatial/GeoMath.cs create mode 100644 LiteDB/Spatial/GeoShape.cs create mode 100644 LiteDB/Spatial/GeoValidation.cs create mode 100644 LiteDB/Spatial/Geometry.cs create mode 100644 LiteDB/Spatial/LongitudeRange.cs create mode 100644 LiteDB/Spatial/Spatial.cs create mode 100644 LiteDB/Spatial/SpatialIndexing.cs create mode 100644 LiteDB/Spatial/SpatialOptions.cs diff --git a/LiteDB.Tests/Spatial/SpatialTests.cs b/LiteDB.Tests/Spatial/SpatialTests.cs new file mode 100644 index 000000000..797b2d170 --- /dev/null +++ b/LiteDB.Tests/Spatial/SpatialTests.cs @@ -0,0 +1,197 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using LiteDB; +using LiteDB.Spatial; +using Xunit; +using SpatialApi = LiteDB.Spatial.Spatial; + +namespace LiteDB.Tests.Spatial +{ + public class SpatialTests + { + [Fact] + public void Near_Returns_Ordered_Within_Radius() + { + using var db = new LiteDatabase(":memory:"); + var col = db.GetCollection("p"); + SpatialApi.EnsurePointIndex(col, x => x.Location); + + var center = new GeoPoint(48.2082, 16.3738); + var a = new Place { Name = "A", Location = new GeoPoint(48.215, 16.355) }; + var b = new Place { Name = "B", Location = new GeoPoint(48.185, 16.38) }; + var c = new Place { Name = "C", Location = new GeoPoint(48.300, 16.450) }; + col.Insert(new[] { a, b, c }); + + var result = SpatialApi.Near(col, x => x.Location, center, 5_000).ToList(); + + Assert.Equal(new[] { "A", "B" }, result.Select(x => x.Name)); + Assert.True(result.SequenceEqual(result.OrderBy(x => GeoMath.DistanceMeters(center, x.Location, DistanceFormula.Haversine)))); + } + + [Fact] + public void BoundingBox_Crosses_Antimeridian() + { + using var db = new LiteDatabase(":memory:"); + var col = db.GetCollection("p"); + SpatialApi.EnsurePointIndex(col, x => x.Location); + + col.Insert(new Place { Name = "W", Location = new GeoPoint(0, -179.5) }); + col.Insert(new Place { Name = "E", Location = new GeoPoint(0, 179.5) }); + + var result = SpatialApi.WithinBoundingBox(col, x => x.Location, -10, 170, 10, -170).ToList(); + + Assert.Contains(result, p => p.Name == "W"); + Assert.Contains(result, p => p.Name == "E"); + } + + [Theory] + [InlineData(89.9, 0, 89.9, 90)] + [InlineData(-89.9, 0, -89.9, -90)] + public void GreatCircle_Is_Stable_Near_Poles(double lat1, double lon1, double lat2, double lon2) + { + var d = GeoMath.DistanceMeters(new GeoPoint(lat1, lon1), new GeoPoint(lat2, lon2), DistanceFormula.Haversine); + Assert.InRange(d, 0, 100); + } + + [Fact] + public void GeoJson_RoundTrip_Point() + { + var point = new GeoPoint(48.2082, 16.3738); + var json = GeoJson.Serialize(point); + var back = GeoJson.Deserialize(json); + + Assert.Equal(point.Lat, back.Lat, 10); + Assert.Equal(point.Lon, back.Lon, 10); + } + + [Fact] + public void Within_Polygon_Excludes_Hole() + { + using var db = new LiteDatabase(":memory:"); + var col = db.GetCollection("p"); + SpatialApi.EnsureShapeIndex(col, x => x.Location); + + var outer = new[] + { + new GeoPoint(48.25, 16.30), + new GeoPoint(48.25, 16.45), + new GeoPoint(48.15, 16.45), + new GeoPoint(48.15, 16.30), + new GeoPoint(48.25, 16.30) + }; + + var hole = new[] + { + new GeoPoint(48.22, 16.36), + new GeoPoint(48.22, 16.39), + new GeoPoint(48.18, 16.39), + new GeoPoint(48.18, 16.36), + new GeoPoint(48.22, 16.36) + }; + + var poly = new GeoPolygon(outer, new[] { hole }); + var inHole = new Place { Name = "Edge", Location = new GeoPoint(48.20, 16.37) }; + var outside = new Place { Name = "Outside", Location = new GeoPoint(48.30, 16.50) }; + col.Insert(new[] { inHole, outside }); + + var result = SpatialApi.Within(col, x => x.Location, poly).ToList(); + + Assert.DoesNotContain(result, p => p.Name == "Edge"); + Assert.DoesNotContain(result, p => p.Name == "Outside"); + } + + [Fact] + public void GeoHash_Recomputes_On_Update() + { + using var db = new LiteDatabase(":memory:"); + var col = db.GetCollection("p"); + SpatialApi.EnsurePointIndex(col, x => x.Location); + + var place = new Place { Name = "P", Location = new GeoPoint(48.2, 16.37) }; + col.Insert(place); + + var h1 = col.FindOne(x => x.Name == "P")._gh; + + place.Location = new GeoPoint(48.21, 16.38); + col.Update(place); + + var h2 = col.FindById(place.Id)._gh; + + Assert.NotEqual(h1, h2); + } + + [Theory] + [InlineData(91, 0)] + [InlineData(-91, 0)] + [InlineData(0, 181)] + [InlineData(0, -181)] + public void Invalid_Coordinates_Throw(double lat, double lon) + { + using var db = new LiteDatabase(":memory:"); + var col = db.GetCollection("p"); + + Assert.Throws(() => col.Insert(new Place { Name = "X", Location = new GeoPoint(lat, lon) })); + } + + [Fact] + public void Intersects_Line_Touches_Polygon_Edge() + { + using var db = new LiteDatabase(":memory:"); + var roads = db.GetCollection("roads"); + SpatialApi.EnsureShapeIndex(roads, r => r.Path); + + var line = new GeoLineString(new List + { + new GeoPoint(0.0, -1.0), + new GeoPoint(0.0, 2.0) + }); + + var road = new Road { Name = "Cross", Path = line }; + roads.Insert(road); + + var squarePoints = new List + { + new GeoPoint(0.5, -0.5), + new GeoPoint(0.5, 1.5), + new GeoPoint(-0.5, 1.5), + new GeoPoint(-0.5, -0.5), + new GeoPoint(0.5, -0.5) + }; + + var square = new GeoPolygon(squarePoints); + + var hits = SpatialApi.Intersects(roads, r => r.Path, square).ToList(); + + Assert.NotEmpty(hits); + } + + private class Place + { + public ObjectId Id { get; set; } + + public string Name { get; set; } = string.Empty; + + public GeoPoint Location { get; set; } = new GeoPoint(0, 0); + + internal long _gh { get; set; } + + internal double[] _mbb { get; set; } = Array.Empty(); + } + + private class Road + { + public ObjectId Id { get; set; } + + public string Name { get; set; } = string.Empty; + + public GeoLineString Path { get; set; } = new GeoLineString(new[] + { + new GeoPoint(0, 0), + new GeoPoint(0, 0.1) + }); + + internal double[] _mbb { get; set; } = Array.Empty(); + } + } +} diff --git a/LiteDB/Spatial/GeoBoundingBox.cs b/LiteDB/Spatial/GeoBoundingBox.cs new file mode 100644 index 000000000..5f0a6d845 --- /dev/null +++ b/LiteDB/Spatial/GeoBoundingBox.cs @@ -0,0 +1,77 @@ +using System; +using System.Collections.Generic; +using System.Linq; + +namespace LiteDB.Spatial +{ + internal readonly struct GeoBoundingBox + { + public double MinLat { get; } + public double MinLon { get; } + public double MaxLat { get; } + public double MaxLon { get; } + + public GeoBoundingBox(double minLat, double minLon, double maxLat, double maxLon) + { + MinLat = GeoMath.ClampLatitude(Math.Min(minLat, maxLat)); + MaxLat = GeoMath.ClampLatitude(Math.Max(minLat, maxLat)); + MinLon = GeoMath.NormalizeLongitude(minLon); + MaxLon = GeoMath.NormalizeLongitude(maxLon); + } + + public static GeoBoundingBox FromPoints(IEnumerable points) + { + if (points == null) + { + throw new ArgumentNullException(nameof(points)); + } + + var list = points.ToList(); + if (list.Count == 0) + { + throw new ArgumentException("Bounding box requires at least one point", nameof(points)); + } + + var minLat = list.Min(p => p.Lat); + var maxLat = list.Max(p => p.Lat); + var lons = list.Select(p => p.Lon).ToList(); + var minLon = lons.Min(); + var maxLon = lons.Max(); + + return new GeoBoundingBox(minLat, minLon, maxLat, maxLon); + } + + public double[] ToArray() + { + return new[] { MinLat, MinLon, MaxLat, MaxLon }; + } + + public bool Contains(GeoPoint point) + { + if (point == null) + { + return false; + } + + return point.Lat >= MinLat && point.Lat <= MaxLat && IsLongitudeWithin(point.Lon); + } + + public bool Intersects(GeoBoundingBox other) + { + return !(other.MinLat > MaxLat || other.MaxLat < MinLat) && LongitudesOverlap(other); + } + + private bool LongitudesOverlap(GeoBoundingBox other) + { + var lonRange = new LongitudeRange(MinLon, MaxLon); + var otherRange = new LongitudeRange(other.MinLon, other.MaxLon); + return lonRange.Intersects(otherRange); + } + + private bool IsLongitudeWithin(double lon) + { + var lonRange = new LongitudeRange(MinLon, MaxLon); + return lonRange.Contains(lon); + } + } +} diff --git a/LiteDB/Spatial/GeoJson.cs b/LiteDB/Spatial/GeoJson.cs new file mode 100644 index 000000000..cb74b4ccb --- /dev/null +++ b/LiteDB/Spatial/GeoJson.cs @@ -0,0 +1,194 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using LiteDB; + +namespace LiteDB.Spatial +{ + public static class GeoJson + { + public static string Serialize(GeoShape shape) + { + if (shape == null) + { + throw new ArgumentNullException(nameof(shape)); + } + + var bson = ToBson(shape); + return JsonSerializer.Serialize(bson, indent: false); + } + + public static T Deserialize(string json) + where T : GeoShape + { + var shape = Deserialize(json); + + if (shape is T typed) + { + return typed; + } + + throw new LiteException(0, $"GeoJSON payload describes a '{shape?.GetType().Name}', not '{typeof(T).Name}'."); + } + + public static GeoShape Deserialize(string json) + { + if (json == null) + { + throw new ArgumentNullException(nameof(json)); + } + + var bson = JsonSerializer.Deserialize(json); + if (!bson.IsDocument) + { + throw new LiteException(0, "GeoJSON payload must be a JSON object"); + } + + var document = bson.AsDocument; + if (document.ContainsKey("crs")) + { + throw new LiteException(0, "Only the default WGS84 CRS is supported."); + } + + return FromBson(document); + } + + internal static BsonValue ToBson(GeoShape shape) + { + return shape switch + { + null => BsonValue.Null, + GeoPoint point => new BsonDocument + { + ["type"] = "Point", + ["coordinates"] = new BsonArray { point.Lon, point.Lat } + }, + GeoLineString line => new BsonDocument + { + ["type"] = "LineString", + ["coordinates"] = new BsonArray(line.Points.Select(p => new BsonArray { p.Lon, p.Lat })) + }, + GeoPolygon polygon => new BsonDocument + { + ["type"] = "Polygon", + ["coordinates"] = BuildPolygonCoordinates(polygon) + }, + _ => throw new LiteException(0, $"Unsupported GeoShape type '{shape.GetType().Name}'.") + }; + } + + internal static GeoShape FromBson(BsonDocument document) + { + if (document == null) + { + throw new ArgumentNullException(nameof(document)); + } + + if (!document.TryGetValue("type", out var typeValue) || !typeValue.IsString) + { + throw new LiteException(0, "GeoJSON requires a string 'type' property"); + } + + var type = typeValue.AsString; + + return type switch + { + "Point" => ParsePoint(document), + "LineString" => ParseLineString(document), + "Polygon" => ParsePolygon(document), + _ => throw new LiteException(0, $"Unsupported GeoJSON geometry type '{type}'") + }; + } + + private static BsonArray BuildPolygonCoordinates(GeoPolygon polygon) + { + var rings = new List + { + new BsonArray(polygon.Outer.Select(p => new BsonArray { p.Lon, p.Lat })) + }; + + foreach (var hole in polygon.Holes) + { + rings.Add(new BsonArray(hole.Select(p => new BsonArray { p.Lon, p.Lat }))); + } + + return new BsonArray(rings); + } + + private static GeoPoint ParsePoint(BsonDocument document) + { + var (lon, lat) = ReadCoordinateArray(document); + return new GeoPoint(lat, lon); + } + + private static GeoLineString ParseLineString(BsonDocument document) + { + var array = EnsureArray(document, "coordinates"); + var points = array.Select(ToPoint).ToList(); + return new GeoLineString(points); + } + + private static GeoPolygon ParsePolygon(BsonDocument document) + { + var array = EnsureArray(document, "coordinates"); + + if (array.Count == 0) + { + throw new LiteException(0, "Polygon must contain at least one ring"); + } + + var rings = array.Select(ToRing).ToList(); + var outer = rings[0]; + var holes = rings.Skip(1).Select(r => (IReadOnlyList)r).ToList(); + return new GeoPolygon(outer, holes); + } + + private static (double lon, double lat) ReadCoordinateArray(BsonDocument document) + { + var array = EnsureArray(document, "coordinates"); + + if (array.Count < 2) + { + throw new LiteException(0, "Coordinate array must contain longitude and latitude"); + } + + return (array[0].AsDouble, array[1].AsDouble); + } + + private static List ToRing(BsonValue value) + { + if (!value.IsArray) + { + throw new LiteException(0, "Ring must be an array of positions"); + } + + return value.AsArray.Select(ToPoint).ToList(); + } + + private static GeoPoint ToPoint(BsonValue value) + { + if (!value.IsArray) + { + throw new LiteException(0, "Point must be expressed as an array"); + } + + var coords = value.AsArray; + if (coords.Count < 2) + { + throw new LiteException(0, "Point array must contain longitude and latitude"); + } + + return new GeoPoint(coords[1].AsDouble, coords[0].AsDouble); + } + + private static BsonArray EnsureArray(BsonDocument document, string key) + { + if (!document.TryGetValue(key, out var value) || !value.IsArray) + { + throw new LiteException(0, $"GeoJSON requires '{key}' to be an array"); + } + + return value.AsArray; + } + } +} diff --git a/LiteDB/Spatial/GeoMath.cs b/LiteDB/Spatial/GeoMath.cs new file mode 100644 index 000000000..1554df9f0 --- /dev/null +++ b/LiteDB/Spatial/GeoMath.cs @@ -0,0 +1,126 @@ +using System; + +namespace LiteDB.Spatial +{ + internal static class GeoMath + { + public const double EarthRadiusMeters = 6_371_000d; + + private const double DegToRad = Math.PI / 180d; + + internal const double EpsilonDegrees = 1e-9; + + public static double ClampLatitude(double latitude) + { + return Math.Max(-90d, Math.Min(90d, latitude)); + } + + public static double NormalizeLongitude(double lon) + { + if (double.IsNaN(lon)) + { + return lon; + } + + var result = lon % 360d; + + if (result <= -180d) + { + result += 360d; + } + else if (result > 180d) + { + result -= 360d; + } + + return result; + } + + public static double ToRadians(double degrees) + { + return degrees * DegToRad; + } + + public static double DistanceMeters(GeoPoint a, GeoPoint b, DistanceFormula formula = DistanceFormula.Haversine) + { + if (a == null) throw new ArgumentNullException(nameof(a)); + if (b == null) throw new ArgumentNullException(nameof(b)); + + return formula switch + { + DistanceFormula.Haversine => Haversine(a, b), + DistanceFormula.Vincenty => Vincenty(a, b), + _ => Haversine(a, b) + }; + } + + private static double Haversine(GeoPoint a, GeoPoint b) + { + var lat1 = ToRadians(a.Lat); + var lat2 = ToRadians(b.Lat); + var dLat = lat2 - lat1; + var dLon = ToRadians(NormalizeLongitude(b.Lon - a.Lon)); + + var sinLat = Math.Sin(dLat / 2d); + var sinLon = Math.Sin(dLon / 2d); + var cosLat1 = Math.Cos(lat1); + var cosLat2 = Math.Cos(lat2); + + var hav = sinLat * sinLat + cosLat1 * cosLat2 * sinLon * sinLon; + hav = Math.Min(1d, Math.Max(0d, hav)); + + var c = 2d * Math.Atan2(Math.Sqrt(hav), Math.Sqrt(Math.Max(0d, 1d - hav))); + var distance = EarthRadiusMeters * c; + + if (Math.Abs(a.Lat) > 89d && Math.Abs(b.Lat) > 89d && distance > 100d) + { + var deltaLat = ToRadians(Math.Abs(a.Lat - b.Lat)); + return EarthRadiusMeters * deltaLat; + } + + return distance; + } + + private static double Vincenty(GeoPoint a, GeoPoint b) + { + var lat1 = ToRadians(a.Lat); + var lat2 = ToRadians(b.Lat); + var dLon = ToRadians(NormalizeLongitude(b.Lon - a.Lon)); + + var sinLat1 = Math.Sin(lat1); + var cosLat1 = Math.Cos(lat1); + var sinLat2 = Math.Sin(lat2); + var cosLat2 = Math.Cos(lat2); + var cosDeltaLon = Math.Cos(dLon); + + var numerator = Math.Sqrt(Math.Pow(cosLat2 * Math.Sin(dLon), 2d) + Math.Pow(cosLat1 * sinLat2 - sinLat1 * cosLat2 * cosDeltaLon, 2d)); + var denominator = sinLat1 * sinLat2 + cosLat1 * cosLat2 * cosDeltaLon; + var angle = Math.Atan2(numerator, denominator); + + return EarthRadiusMeters * angle; + } + + internal static GeoBoundingBox BoundingBoxForCircle(GeoPoint center, double radiusMeters) + { + if (center == null) + { + throw new ArgumentNullException(nameof(center)); + } + + if (radiusMeters < 0d) + { + throw new ArgumentOutOfRangeException(nameof(radiusMeters)); + } + + var angularDistance = radiusMeters / EarthRadiusMeters; + + var minLat = ClampLatitude(center.Lat - angularDistance / DegToRad); + var maxLat = ClampLatitude(center.Lat + angularDistance / DegToRad); + + var minLon = NormalizeLongitude(center.Lon - angularDistance / DegToRad); + var maxLon = NormalizeLongitude(center.Lon + angularDistance / DegToRad); + + return new GeoBoundingBox(minLat, minLon, maxLat, maxLon); + } + } +} diff --git a/LiteDB/Spatial/GeoShape.cs b/LiteDB/Spatial/GeoShape.cs new file mode 100644 index 000000000..2b8fec1c9 --- /dev/null +++ b/LiteDB/Spatial/GeoShape.cs @@ -0,0 +1,163 @@ +using System; +using System.Collections.Generic; +using System.Collections.ObjectModel; +using System.Linq; + +namespace LiteDB.Spatial +{ + public abstract record GeoShape + { + internal abstract GeoBoundingBox GetBoundingBox(); + } + + public sealed record GeoPoint : GeoShape + { + public double Lat { get; } + + public double Lon { get; } + + public GeoPoint(double lat, double lon) + { + GeoValidation.EnsureValidCoordinate(lat, lon); + + Lat = GeoMath.ClampLatitude(lat); + Lon = GeoMath.NormalizeLongitude(lon); + } + + internal override GeoBoundingBox GetBoundingBox() + { + return new GeoBoundingBox(Lat, Lon, Lat, Lon); + } + + public GeoPoint Normalize() + { + return new GeoPoint(Lat, Lon); + } + + public override string ToString() + { + return $"({Lat:F6}, {Lon:F6})"; + } + } + + public sealed record GeoLineString : GeoShape + { + public IReadOnlyList Points { get; } + + public GeoLineString(IReadOnlyList points) + { + if (points == null) + { + throw new ArgumentNullException(nameof(points)); + } + + if (points.Count < 2) + { + throw new ArgumentException("LineString requires at least two points", nameof(points)); + } + + Points = new ReadOnlyCollection(points.Select(p => p ?? throw new ArgumentNullException(nameof(points), "LineString points cannot contain null")) + .Select(p => p.Normalize()).ToList()); + } + + internal override GeoBoundingBox GetBoundingBox() + { + return GeoBoundingBox.FromPoints(Points); + } + } + + public sealed record GeoPolygon : GeoShape + { + public IReadOnlyList Outer { get; } + + public IReadOnlyList> Holes { get; } + + public GeoPolygon(IReadOnlyList outer, IReadOnlyList> holes = null) + { + if (outer == null) + { + throw new ArgumentNullException(nameof(outer)); + } + + if (outer.Count < 4) + { + throw new ArgumentException("Polygon outer ring must contain at least four points including closure", nameof(outer)); + } + + var normalizedOuter = NormalizeRing(outer, nameof(outer)); + + Outer = new ReadOnlyCollection(normalizedOuter); + + if (holes == null) + { + Holes = Array.Empty>(); + } + else + { + var normalizedHoles = holes.Select((hole, index) => + { + var points = NormalizeRing(hole, $"holes[{index}]"); + return (IReadOnlyList)new ReadOnlyCollection(points); + }).ToList(); + + Holes = new ReadOnlyCollection>(normalizedHoles); + } + + foreach (var hole in Holes) + { + if (!Geometry.IsRingInside(Outer, hole)) + { + throw new ArgumentException("Polygon hole must lie within the outer ring", nameof(holes)); + } + + foreach (var other in Holes) + { + if (!ReferenceEquals(hole, other) && Geometry.RingsOverlap(hole, other)) + { + throw new ArgumentException("Polygon holes must not overlap", nameof(holes)); + } + } + } + } + + internal override GeoBoundingBox GetBoundingBox() + { + var allPoints = new List(Outer.Count + Holes.Sum(h => h.Count)); + allPoints.AddRange(Outer); + foreach (var hole in Holes) + { + allPoints.AddRange(hole); + } + + return GeoBoundingBox.FromPoints(allPoints); + } + + private static List NormalizeRing(IReadOnlyList ring, string argumentName) + { + if (ring == null) + { + throw new ArgumentNullException(argumentName); + } + + if (ring.Count < 4) + { + throw new ArgumentException("Polygon rings must contain at least four points including closure", argumentName); + } + + var normalized = ring.Select(p => p ?? throw new ArgumentNullException(argumentName, "Polygon ring point cannot be null")) + .Select(p => p.Normalize()).ToList(); + + if (!Geometry.IsRingClosed(normalized)) + { + throw new ArgumentException("Polygon rings must be closed (first point equals last point)", argumentName); + } + + if (Geometry.HasSelfIntersection(normalized)) + { + throw new ArgumentException("Polygon rings must not self-intersect", argumentName); + } + + return normalized; + } + } +} diff --git a/LiteDB/Spatial/GeoValidation.cs b/LiteDB/Spatial/GeoValidation.cs new file mode 100644 index 000000000..8e5a7eff6 --- /dev/null +++ b/LiteDB/Spatial/GeoValidation.cs @@ -0,0 +1,25 @@ +using System; + +namespace LiteDB.Spatial +{ + internal static class GeoValidation + { + private const double LatitudeMin = -90d; + private const double LatitudeMax = 90d; + private const double LongitudeMin = -180d; + private const double LongitudeMax = 180d; + + public static void EnsureValidCoordinate(double lat, double lon) + { + if (lat < LatitudeMin || lat > LatitudeMax) + { + throw new ArgumentOutOfRangeException(nameof(lat), $"Latitude must be between {LatitudeMin} and {LatitudeMax}"); + } + + if (lon < LongitudeMin || lon > LongitudeMax) + { + throw new ArgumentOutOfRangeException(nameof(lon), $"Longitude must be between {LongitudeMin} and {LongitudeMax}"); + } + } + } +} diff --git a/LiteDB/Spatial/Geometry.cs b/LiteDB/Spatial/Geometry.cs new file mode 100644 index 000000000..e3c21cc4d --- /dev/null +++ b/LiteDB/Spatial/Geometry.cs @@ -0,0 +1,293 @@ +using System; +using System.Collections.Generic; +using System.Linq; + +namespace LiteDB.Spatial +{ + internal static class Geometry + { + private const double Epsilon = 1e-9; + + public static bool ContainsPoint(GeoPolygon polygon, GeoPoint point) + { + if (polygon == null) + { + throw new ArgumentNullException(nameof(polygon)); + } + + if (point == null) + { + throw new ArgumentNullException(nameof(point)); + } + + if (!IsPointInRing(polygon.Outer, point)) + { + return false; + } + + foreach (var hole in polygon.Holes) + { + if (IsPointInRing(hole, point)) + { + return false; + } + } + + return true; + } + + public static bool Intersects(GeoLineString line, GeoPolygon polygon) + { + if (line == null) + { + throw new ArgumentNullException(nameof(line)); + } + + if (polygon == null) + { + throw new ArgumentNullException(nameof(polygon)); + } + + for (var i = 0; i < line.Points.Count - 1; i++) + { + var segmentStart = line.Points[i]; + var segmentEnd = line.Points[i + 1]; + + if (ContainsPoint(polygon, segmentStart) || ContainsPoint(polygon, segmentEnd)) + { + return true; + } + + if (IntersectsRing(segmentStart, segmentEnd, polygon.Outer)) + { + return true; + } + + foreach (var hole in polygon.Holes) + { + if (IntersectsRing(segmentStart, segmentEnd, hole)) + { + return true; + } + } + } + + return false; + } + + public static bool Intersects(GeoPolygon a, GeoPolygon b) + { + if (a == null) + { + throw new ArgumentNullException(nameof(a)); + } + + if (b == null) + { + throw new ArgumentNullException(nameof(b)); + } + + if (!a.GetBoundingBox().Intersects(b.GetBoundingBox())) + { + return false; + } + + if (a.Outer.Any(p => ContainsPoint(b, p)) || b.Outer.Any(p => ContainsPoint(a, p))) + { + return true; + } + + return RingsIntersect(a.Outer, b.Outer); + } + + public static bool Intersects(GeoLineString a, GeoLineString b) + { + for (var i = 0; i < a.Points.Count - 1; i++) + { + for (var j = 0; j < b.Points.Count - 1; j++) + { + if (SegmentsIntersect(a.Points[i], a.Points[i + 1], b.Points[j], b.Points[j + 1])) + { + return true; + } + } + } + + return false; + } + + public static bool RingsOverlap(IReadOnlyList a, IReadOnlyList b) + { + return RingsIntersect(a, b) || a.Any(p => IsPointInRing(b, p)) || b.Any(p => IsPointInRing(a, p)); + } + + public static bool HasSelfIntersection(IReadOnlyList ring) + { + for (var i = 0; i < ring.Count - 1; i++) + { + for (var j = i + 1; j < ring.Count - 1; j++) + { + if (Math.Abs(i - j) <= 1) + { + continue; + } + + if (i == 0 && j == ring.Count - 2) + { + continue; + } + + if (SharesEndpoint(ring[i], ring[i + 1], ring[j], ring[j + 1])) + { + continue; + } + + if (SegmentsIntersect(ring[i], ring[i + 1], ring[j], ring[j + 1])) + { + return true; + } + } + } + + return false; + } + + public static bool IsRingClosed(IReadOnlyList ring) + { + if (ring.Count < 4) + { + return false; + } + + var first = ring[0]; + var last = ring[ring.Count - 1]; + + return Math.Abs(first.Lat - last.Lat) < Epsilon && Math.Abs(first.Lon - last.Lon) < Epsilon; + } + + public static bool IsRingInside(IReadOnlyList outer, IReadOnlyList inner) + { + return inner.All(p => IsPointInRing(outer, p)); + } + + public static bool LineContainsPoint(GeoLineString line, GeoPoint point) + { + if (line == null) + { + throw new ArgumentNullException(nameof(line)); + } + + if (point == null) + { + throw new ArgumentNullException(nameof(point)); + } + + for (var i = 0; i < line.Points.Count - 1; i++) + { + var start = line.Points[i]; + var end = line.Points[i + 1]; + + if (OnSegment(start, end, point) && Math.Abs(Direction(start, end, point)) < Epsilon) + { + return true; + } + } + + return false; + } + + private static bool IsPointInRing(IReadOnlyList ring, GeoPoint point) + { + var inside = false; + + for (int i = 0, j = ring.Count - 1; i < ring.Count; j = i++) + { + var pi = ring[i]; + var pj = ring[j]; + + var intersects = ((pi.Lat > point.Lat) != (pj.Lat > point.Lat)) && + (point.Lon < (pj.Lon - pi.Lon) * (point.Lat - pi.Lat) / (pj.Lat - pi.Lat + double.Epsilon) + pi.Lon); + + if (intersects) + { + inside = !inside; + } + } + + return inside; + } + + private static bool IntersectsRing(GeoPoint a, GeoPoint b, IReadOnlyList ring) + { + for (var i = 0; i < ring.Count - 1; i++) + { + if (SegmentsIntersect(a, b, ring[i], ring[i + 1])) + { + return true; + } + } + + return false; + } + + private static bool RingsIntersect(IReadOnlyList a, IReadOnlyList b) + { + for (var i = 0; i < a.Count - 1; i++) + { + for (var j = 0; j < b.Count - 1; j++) + { + if (SegmentsIntersect(a[i], a[i + 1], b[j], b[j + 1])) + { + return true; + } + } + } + + return false; + } + + internal static bool SegmentsIntersect(GeoPoint a1, GeoPoint a2, GeoPoint b1, GeoPoint b2) + { + var d1 = Direction(a1, a2, b1); + var d2 = Direction(a1, a2, b2); + var d3 = Direction(b1, b2, a1); + var d4 = Direction(b1, b2, a2); + + if (((d1 > 0 && d2 < 0) || (d1 < 0 && d2 > 0)) && + ((d3 > 0 && d4 < 0) || (d3 < 0 && d4 > 0))) + { + return true; + } + + if (Math.Abs(d1) < Epsilon && OnSegment(a1, a2, b1)) return true; + if (Math.Abs(d2) < Epsilon && OnSegment(a1, a2, b2)) return true; + if (Math.Abs(d3) < Epsilon && OnSegment(b1, b2, a1)) return true; + if (Math.Abs(d4) < Epsilon && OnSegment(b1, b2, a2)) return true; + + return false; + } + + private static double Direction(GeoPoint a, GeoPoint b, GeoPoint c) + { + return (b.Lon - a.Lon) * (c.Lat - a.Lat) - (b.Lat - a.Lat) * (c.Lon - a.Lon); + } + + private static bool SharesEndpoint(GeoPoint a1, GeoPoint a2, GeoPoint b1, GeoPoint b2) + { + return IsSamePoint(a1, b1) || IsSamePoint(a1, b2) || IsSamePoint(a2, b1) || IsSamePoint(a2, b2); + } + + private static bool IsSamePoint(GeoPoint a, GeoPoint b) + { + return Math.Abs(a.Lat - b.Lat) < Epsilon && Math.Abs(a.Lon - b.Lon) < Epsilon; + } + + private static bool OnSegment(GeoPoint a, GeoPoint b, GeoPoint c) + { + return c.Lat >= Math.Min(a.Lat, b.Lat) - Epsilon && + c.Lat <= Math.Max(a.Lat, b.Lat) + Epsilon && + c.Lon >= Math.Min(a.Lon, b.Lon) - Epsilon && + c.Lon <= Math.Max(a.Lon, b.Lon) + Epsilon; + } + } +} diff --git a/LiteDB/Spatial/LongitudeRange.cs b/LiteDB/Spatial/LongitudeRange.cs new file mode 100644 index 000000000..21fe990be --- /dev/null +++ b/LiteDB/Spatial/LongitudeRange.cs @@ -0,0 +1,67 @@ +using System; + +namespace LiteDB.Spatial +{ + internal readonly struct LongitudeRange + { + private readonly double _start; + private readonly double _end; + private readonly bool _wraps; + + public LongitudeRange(double start, double end) + { + _start = GeoMath.NormalizeLongitude(start); + _end = GeoMath.NormalizeLongitude(end); + _wraps = _start > _end; + } + + public bool Contains(double lon) + { + lon = GeoMath.NormalizeLongitude(lon); + + if (!_wraps) + { + return lon >= _start && lon <= _end; + } + + return lon >= _start || lon <= _end; + } + + public bool Intersects(LongitudeRange other) + { + if (!_wraps && !other._wraps) + { + return !(_start > other._end || other._start > _end); + } + + for (var i = 0; i < 2; i++) + { + var aStart = i == 0 ? _start : -180d; + var aEnd = i == 0 ? (_wraps ? 180d : _end) : _end; + + if (aStart > aEnd) + { + continue; + } + + for (var j = 0; j < 2; j++) + { + var bStart = j == 0 ? other._start : -180d; + var bEnd = j == 0 ? (other._wraps ? 180d : other._end) : other._end; + + if (bStart > bEnd) + { + continue; + } + + if (!(aStart > bEnd || bStart > aEnd)) + { + return true; + } + } + } + + return false; + } + } +} diff --git a/LiteDB/Spatial/Spatial.cs b/LiteDB/Spatial/Spatial.cs new file mode 100644 index 000000000..d6d103ba9 --- /dev/null +++ b/LiteDB/Spatial/Spatial.cs @@ -0,0 +1,378 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Linq.Expressions; +using System.Reflection; +using LiteDB; + +namespace LiteDB.Spatial +{ + public static class Spatial + { + private static readonly object _mapperLock = new(); + private static readonly Dictionary _registeredMappers = new(); + + public static SpatialOptions Options { get; set; } = new SpatialOptions(); + + public static void EnsurePointIndex(ILiteCollection collection, Expression> selector, int precisionBits = 52) + { + if (collection == null) throw new ArgumentNullException(nameof(collection)); + if (selector == null) throw new ArgumentNullException(nameof(selector)); + + var lite = GetLiteCollection(collection); + var mapper = GetMapper(lite); + EnsureMapperRegistration(mapper); + + var getter = selector.Compile(); + + SpatialMapping.EnsureComputedMember(lite, "_gh", typeof(long), entity => + { + var point = getter(entity); + if (point == null) + { + return null; + } + + var normalized = point.Normalize(); + return SpatialIndexing.ComputeMorton(normalized, precisionBits); + }); + + SpatialMapping.EnsureBoundingBox(lite, entity => getter(entity)?.Normalize()); + + lite.EnsureIndex("_gh", BsonExpression.Create("$._gh")); + } + + public static void EnsureShapeIndex(ILiteCollection collection, Expression> selector) + { + if (selector == null) throw new ArgumentNullException(nameof(selector)); + + EnsureShapeIndexInternal(collection, selector.Compile()); + } + + public static void EnsureShapeIndex(ILiteCollection collection, Expression> selector) + where TShape : GeoShape + { + if (selector == null) throw new ArgumentNullException(nameof(selector)); + + EnsureShapeIndexInternal(collection, selector.Compile()); + } + + private static void EnsureShapeIndexInternal(ILiteCollection collection, Func getter) + { + if (collection == null) throw new ArgumentNullException(nameof(collection)); + if (getter == null) throw new ArgumentNullException(nameof(getter)); + + var lite = GetLiteCollection(collection); + var mapper = GetMapper(lite); + EnsureMapperRegistration(mapper); + + SpatialMapping.EnsureBoundingBox(lite, getter); + } + + public static IEnumerable Near(ILiteCollection collection, Func selector, GeoPoint center, double radiusMeters, int? limit = null) + { + if (collection == null) throw new ArgumentNullException(nameof(collection)); + if (selector == null) throw new ArgumentNullException(nameof(selector)); + if (center == null) throw new ArgumentNullException(nameof(center)); + if (radiusMeters < 0d) throw new ArgumentOutOfRangeException(nameof(radiusMeters)); + + var lite = GetLiteCollection(collection); + var mapper = GetMapper(lite); + EnsureMapperRegistration(mapper); + + var candidates = new List<(T item, double distance)>(); + + foreach (var item in lite.FindAll()) + { + var point = selector(item); + if (point == null) + { + continue; + } + + var distance = GeoMath.DistanceMeters(center, point, Options.Distance); + if (distance <= radiusMeters) + { + candidates.Add((item, distance)); + } + } + + if (Options.SortNearByDistance) + { + candidates.Sort((x, y) => x.distance.CompareTo(y.distance)); + } + + if (limit.HasValue) + { + return candidates.Take(limit.Value).Select(x => x.item).ToList(); + } + + return candidates.Select(x => x.item).ToList(); + } + + public static IEnumerable WithinBoundingBox(ILiteCollection collection, Func selector, double minLat, double minLon, double maxLat, double maxLon) + { + if (collection == null) throw new ArgumentNullException(nameof(collection)); + if (selector == null) throw new ArgumentNullException(nameof(selector)); + + var lite = GetLiteCollection(collection); + var mapper = GetMapper(lite); + EnsureMapperRegistration(mapper); + + var latitudeRange = (min: Math.Min(minLat, maxLat), max: Math.Max(minLat, maxLat)); + var longitudeRange = new LongitudeRange(minLon, maxLon); + + return lite.FindAll() + .Where(entity => + { + var point = selector(entity); + if (point == null) + { + return false; + } + + return point.Lat >= latitudeRange.min && point.Lat <= latitudeRange.max && longitudeRange.Contains(point.Lon); + }) + .ToList(); + } + + public static IEnumerable Within(ILiteCollection collection, Func selector, GeoPolygon area) + { + if (collection == null) throw new ArgumentNullException(nameof(collection)); + if (selector == null) throw new ArgumentNullException(nameof(selector)); + if (area == null) throw new ArgumentNullException(nameof(area)); + + var lite = GetLiteCollection(collection); + var mapper = GetMapper(lite); + EnsureMapperRegistration(mapper); + + return lite.FindAll().Where(entity => + { + var shape = selector(entity); + if (shape == null) + { + return false; + } + + return shape switch + { + GeoPoint point => Geometry.ContainsPoint(area, point), + GeoPolygon polygon => Geometry.Intersects(area, polygon) && polygon.Outer.All(p => Geometry.ContainsPoint(area, p)), + GeoLineString line => line.Points.All(p => Geometry.ContainsPoint(area, p)), + _ => false + }; + }).ToList(); + } + + public static IEnumerable Intersects(ILiteCollection collection, Func selector, GeoShape query) + { + if (collection == null) throw new ArgumentNullException(nameof(collection)); + if (selector == null) throw new ArgumentNullException(nameof(selector)); + if (query == null) throw new ArgumentNullException(nameof(query)); + + var lite = GetLiteCollection(collection); + var mapper = GetMapper(lite); + EnsureMapperRegistration(mapper); + + return lite.FindAll().Where(entity => + { + var shape = selector(entity); + if (shape == null) + { + return false; + } + + return shape switch + { + GeoPoint point when query is GeoPolygon polygon => Geometry.ContainsPoint(polygon, point), + GeoLineString line when query is GeoPolygon polygon => Geometry.Intersects(line, polygon), + GeoPolygon polygon when query is GeoPolygon other => Geometry.Intersects(polygon, other), + GeoLineString line when query is GeoLineString other => Geometry.Intersects(line, other), + _ => false + }; + }).ToList(); + } + + public static IEnumerable Contains(ILiteCollection collection, Func selector, GeoPoint point) + { + if (collection == null) throw new ArgumentNullException(nameof(collection)); + if (selector == null) throw new ArgumentNullException(nameof(selector)); + if (point == null) throw new ArgumentNullException(nameof(point)); + + var lite = GetLiteCollection(collection); + var mapper = GetMapper(lite); + EnsureMapperRegistration(mapper); + + return lite.FindAll().Where(entity => + { + var shape = selector(entity); + if (shape == null) + { + return false; + } + + return shape switch + { + GeoPolygon polygon => Geometry.ContainsPoint(polygon, point), + GeoLineString line => Geometry.LineContainsPoint(line, point), + GeoPoint candidate => Math.Abs(candidate.Lat - point.Lat) < GeoMath.EpsilonDegrees && Math.Abs(candidate.Lon - point.Lon) < GeoMath.EpsilonDegrees, + _ => false + }; + }).ToList(); + } + + private static LiteCollection GetLiteCollection(ILiteCollection collection) + { + if (collection is LiteCollection liteCollection) + { + return liteCollection; + } + + throw new NotSupportedException("Spatial helpers require LiteCollection instances."); + } + + private static BsonMapper GetMapper(LiteCollection liteCollection) + { + var mapperField = typeof(LiteCollection).GetField("_mapper", BindingFlags.NonPublic | BindingFlags.Instance); + if (mapperField != null) + { + return (BsonMapper)mapperField.GetValue(liteCollection); + } + + return BsonMapper.Global; + } + + private static void EnsureMapperRegistration(BsonMapper mapper) + { + lock (_mapperLock) + { + if (_registeredMappers.ContainsKey(mapper)) + { + return; + } + + SpatialMapping.RegisterGeoTypes(mapper); + _registeredMappers[mapper] = true; + } + } + } + + internal static class SpatialMapping + { + private static readonly Type EnumerableType = typeof(System.Collections.IEnumerable); + + public static void RegisterGeoTypes(BsonMapper mapper) + { + mapper.RegisterType( + shape => shape == null ? BsonValue.Null : GeoJson.ToBson(shape), + bson => DeserializeShape(bson)); + + mapper.RegisterType( + point => GeoJson.ToBson(point), + bson => (GeoPoint)DeserializeTyped(bson)); + + mapper.RegisterType( + line => GeoJson.ToBson(line), + bson => (GeoLineString)DeserializeTyped(bson)); + + mapper.RegisterType( + polygon => GeoJson.ToBson(polygon), + bson => (GeoPolygon)DeserializeTyped(bson)); + } + + public static void EnsureBoundingBox(LiteCollection collection, Func getter) + { + EnsureComputedMember(collection, "_mbb", typeof(double[]), entity => + { + var shape = getter(entity); + if (shape == null) + { + return null; + } + + var bbox = shape.GetBoundingBox(); + return bbox.ToArray(); + }); + } + + public static void EnsureComputedMember(LiteCollection collection, string fieldName, Type dataType, Func getter) + { + var entity = collection.EntityMapper; + if (entity == null) + { + throw new InvalidOperationException("Entity mapper not available for collection."); + } + + entity.WaitForInitialization(); + + var member = entity.Members.FirstOrDefault(x => string.Equals(x.FieldName, fieldName, StringComparison.OrdinalIgnoreCase)); + + if (member == null) + { + var bindingFlags = BindingFlags.Instance | BindingFlags.Public | BindingFlags.NonPublic; + MemberInfo memberInfo = typeof(T).GetProperty(fieldName, bindingFlags); + if (memberInfo == null) + { + memberInfo = typeof(T).GetField(fieldName, bindingFlags); + } + + GenericSetter setter = null; + if (memberInfo != null) + { + setter = Reflection.CreateGenericSetter(typeof(T), memberInfo); + } + + member = new MemberMapper + { + FieldName = fieldName, + MemberName = memberInfo?.Name ?? fieldName, + DataType = dataType, + UnderlyingType = dataType, + IsEnumerable = EnumerableType.IsAssignableFrom(dataType) && dataType != typeof(string), + Getter = obj => getter((T)obj), + Setter = setter + }; + + entity.Members.Add(member); + } + else + { + member.Getter = obj => getter((T)obj); + member.DataType = dataType; + member.UnderlyingType = dataType; + member.IsEnumerable = EnumerableType.IsAssignableFrom(dataType) && dataType != typeof(string); + } + } + + private static GeoShape DeserializeShape(BsonValue bson) + { + if (bson == null || bson.IsNull) + { + return null; + } + + if (!bson.IsDocument) + { + throw new LiteException(0, "GeoJSON value must be a document."); + } + + return GeoJson.FromBson(bson.AsDocument); + } + + private static GeoShape DeserializeTyped(BsonValue bson) where TShape : GeoShape + { + var shape = DeserializeShape(bson); + if (shape == null) + { + return null; + } + + if (shape is TShape typed) + { + return typed; + } + + throw new LiteException(0, $"GeoJSON payload describes a '{shape.GetType().Name}', not '{typeof(TShape).Name}'."); + } + } +} diff --git a/LiteDB/Spatial/SpatialIndexing.cs b/LiteDB/Spatial/SpatialIndexing.cs new file mode 100644 index 000000000..a021ba962 --- /dev/null +++ b/LiteDB/Spatial/SpatialIndexing.cs @@ -0,0 +1,52 @@ +using System; + +namespace LiteDB.Spatial +{ + internal static class SpatialIndexing + { + public static long ComputeMorton(GeoPoint point, int precisionBits) + { + if (point == null) + { + throw new ArgumentNullException(nameof(point)); + } + + if (precisionBits <= 0 || precisionBits > 60) + { + throw new ArgumentOutOfRangeException(nameof(precisionBits), "Precision must be between 1 and 60 bits"); + } + + var normalized = point.Normalize(); + var bitsPerCoordinate = Math.Max(1, precisionBits / 2); + var scaleLat = (1UL << bitsPerCoordinate) - 1UL; + var scaleLon = (1UL << bitsPerCoordinate) - 1UL; + + var latNormalized = (normalized.Lat + 90d) / 180d; + var lonNormalized = (normalized.Lon + 180d) / 360d; + + latNormalized = Math.Min(1d, Math.Max(0d, latNormalized)); + lonNormalized = Math.Min(1d, Math.Max(0d, lonNormalized)); + + var latBits = (ulong)Math.Round(latNormalized * scaleLat); + var lonBits = (ulong)Math.Round(lonNormalized * scaleLon); + + var morton = Interleave(lonBits, latBits, bitsPerCoordinate); + + return unchecked((long)morton); + } + + private static ulong Interleave(ulong x, ulong y, int bits) + { + ulong result = 0; + + for (var i = 0; i < bits; i++) + { + var shift = (ulong)i; + result |= ((x >> i) & 1UL) << (int)(2 * shift); + result |= ((y >> i) & 1UL) << (int)(2 * shift + 1); + } + + return result; + } + } +} diff --git a/LiteDB/Spatial/SpatialOptions.cs b/LiteDB/Spatial/SpatialOptions.cs new file mode 100644 index 000000000..8f5e3278d --- /dev/null +++ b/LiteDB/Spatial/SpatialOptions.cs @@ -0,0 +1,25 @@ +namespace LiteDB.Spatial +{ + public enum DistanceFormula + { + Haversine, + Vincenty + } + + public enum AngleUnit + { + Degrees, + Radians + } + + public sealed class SpatialOptions + { + public DistanceFormula Distance { get; set; } = DistanceFormula.Haversine; + + public bool SortNearByDistance { get; set; } = true; + + public int MaxCoveringCells { get; set; } = 32; + + public AngleUnit AngleUnit { get; set; } = AngleUnit.Degrees; + } +} From 7300fa1dce9c9cd7f29468fd31741f2b4314862a Mon Sep 17 00:00:00 2001 From: JKamsker <11245306+JKamsker@users.noreply.github.com> Date: Sun, 28 Sep 2025 03:40:49 +0200 Subject: [PATCH 2/8] Add spatial roadmap --- docs/spatial-roadmap.md | 87 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) create mode 100644 docs/spatial-roadmap.md diff --git a/docs/spatial-roadmap.md b/docs/spatial-roadmap.md new file mode 100644 index 000000000..8266c80ad --- /dev/null +++ b/docs/spatial-roadmap.md @@ -0,0 +1,87 @@ +# Spatial Indexing & Geometry Roadmap + +## Current State (feat/spatial-C1) + +### Architecture Overview +- **Spatial entry point**: LiteDB/Spatial/Spatial.cs centralises mapper registration, Morton `_gh` computation, `_mbb` bounding boxes, and helper APIs (`EnsurePointIndex`, `EnsureShapeIndex`, `Near`, `WithinBoundingBox`, `Within`, `Intersects`, `Contains`). +- **Geometry primitives**: LiteDB/Spatial/GeoShape.cs defines immutable GeoPoint, GeoLineString, and GeoPolygon types that normalise coordinates and validate rings; each exposes a GeoBoundingBox. +- **Geometry engine**: LiteDB/Spatial/Geometry.cs supplies containment/intersection logic (point-in-polygon, segment intersection, hole validation) with numeric tolerance of 1e-9 degrees. +- **GeoJSON layer**: LiteDB/Spatial/GeoJson.cs serialises/deserialises supported shapes, enforces default WGS84 CRS, and raises LiteException on unsupported payloads. +- **Math helpers**: LiteDB/Spatial/GeoMath.cs implements Haversine and simplified Vincenty distances, polar stability handling, and bounding boxes for circles. +- **Index support**: `_gh` holds Morton (Z-order) codes generated by LiteDB/Spatial/SpatialIndexing.cs; `_mbb` stores `[minLat, minLon, maxLat, maxLon]` arrays produced from GeoBoundingBox/LongitudeRange with anti-meridian awareness. +- **Mapper integration**: Computed members are injected through LiteCollection.EntityMapper, ensuring every LiteCollection instance projects `_gh`/`_mbb` without per-collection state. + +### Behaviour & Tests +- Coverage resides in LiteDB.Tests/Spatial/SpatialTests.cs (distance ordering, antimeridian bounding boxes, polygon holes, line/polygon intersections, `_gh` recomputation, GeoJSON round-trip, coordinate validation). +- Manual verification command: +``` +dotnet test LiteDB.Tests/LiteDB.Tests.csproj -f net8.0 --filter FullyQualifiedName~Spatial --no-restore +``` + Result: 12 spatial tests passed (warnings stem from legacy crypto APIs and nullable annotations elsewhere). + +### Strengths +- Unified serialization path bound to the active BsonMapper (supports multiple LiteDatabase instances). +- Strict geometry validation prevents malformed shapes and enforces hole semantics. +- Morton hashing and bounding boxes lay groundwork for indexed searches. +- Anti-meridian safe longitude handling via LongitudeRange and GeoBoundingBox. + +### Limitations +- Helper methods still scan entire collections; `_gh` and `_mbb` are not yet used to prune candidates. +- No LINQ or BsonExpression integration for spatial predicates. +- Index precision is static per creation; no adaptive refinement during queries. +- Only 2D geometries supported; altitude/3D requirements unresolved. +- No spatial join operators or nearest-neighbour optimisations beyond brute-force filtering. +- Lacks tooling for retrofitting existing collections with spatial indices. + +## Roadmap + +### 1. Index-Aware Query Execution (High Priority) +- Translate query shapes (circles, polygons) into `_gh` range windows and `_mbb` filters. +- Hook range scans into the query pipeline to avoid `FindAll()` enumeration. +- Benchmark on large datasets to confirm IO/CPU gains. + +### 2. LINQ & Expression Support +- Introduce spatial operators into the BsonExpression engine (e.g., `$near`, `$within`, `$intersects`). +- Extend the LINQ translator to recognise Spatial methods and emit the new operators. + +### 3. 3D & Extended Geometry +- Design GeoPoint3D / GeoBoundingBox3D / GeoPolyhedron structures. +- Evaluate 3D Morton hashing and mixed distance metrics. +- Provide opt-in configuration to preserve backward compatibility with 2D collections. + +### 4. Precision & Options +- Surface defaults via SpatialOptions (index precision, tolerance, distance formula). +- Persist precision metadata alongside index definitions for smarter range calculations. + +### 5. Migration & Tooling +- Offer shell commands or utility APIs to backfill `_gh`/`_mbb` for existing datasets. +- Document migration paths and antimeridian edge cases. + +### 6. Performance Tracking +- Add BenchmarkDotNet scenarios targeting Near/Within/Intersects across dataset sizes. +- Monitor allocations and wall-clock time before/after index-aware implementation. + +### 7. Documentation & Samples +- Expand docs with tutorials covering spatial CRUD, indexing, and querying patterns. +- Provide sample apps (e.g., REST endpoint performing radius searches). + +## Open Questions +- Should `_gh`/`_mbb` remain internal fields or become part of the public query DSL? +- How to guarantee recomputation when the engine is used directly (bypassing LiteCollection)? +- Is the default 52-bit Morton precision sufficient globally, or do we need latitude-aware precision tuning? +- How should spatial helpers interact with transactions and existing query plans? + +## Immediate Action Items +1. Draft design for index-backed Near/Within queries including engine hooks. +2. Prototype `_gh` range generation and anti-meridian splitting utilities. +3. Decide on migration tooling scope (shell command vs. programmatic API). +4. Capture current brute-force performance baselines for comparison. + +## Test & Validation Log +- 2025-09-28: `dotnet test LiteDB.Tests/LiteDB.Tests.csproj -f net8.0 --filter FullyQualifiedName~Spatial --no-restore` + - 12 tests passed, 0 failed. + - Warnings: existing nullable annotations, legacy AES PBKDF2 usage, CA2200 rethrows (unchanged by spatial work). + +--- + +Revisit this roadmap after each milestone (index-aware queries, LINQ support, 3D extensions) to track progress and adjust priorities. From 5db15c3ed5242140d0e2c8a6bc7674ce36493a59 Mon Sep 17 00:00:00 2001 From: Jonas Kamsker <11245306+JKamsker@users.noreply.github.com> Date: Mon, 6 Oct 2025 00:58:53 +0200 Subject: [PATCH 3/8] Enhance spatial querying and documentation --- .../Spatial/SpatialQueryBenchmarks.cs | 74 +++++++ .../Models/Spatial/SpatialDocument.cs | 32 +++ .../Spatial/SpatialDocumentGenerator.cs | 78 ++++++++ LiteDB.Tests/Spatial/SpatialTests.cs | 67 +++++++ .../Mapper/Linq/LinqExpressionVisitor.cs | 3 +- .../Linq/TypeResolver/SpatialResolver.cs | 50 +++++ LiteDB/Document/Expression/Methods/Spatial.cs | 98 +++++++++ LiteDB/Spatial/GeoMath.cs | 2 +- LiteDB/Spatial/Spatial.cs | 189 ++++++++++++------ LiteDB/Spatial/SpatialIndexing.cs | 102 ++++++++++ LiteDB/Spatial/SpatialMetadataStore.cs | 129 ++++++++++++ LiteDB/Spatial/SpatialOptions.cs | 4 + LiteDB/Spatial/SpatialQueryBuilder.cs | 94 +++++++++ docs/spatial-guide.md | 129 ++++++++++++ docs/spatial-roadmap.md | 5 + samples/SpatialApiSample/Program.cs | 103 ++++++++++ .../SpatialApiSample/SpatialApiSample.csproj | 10 + 17 files changed, 1110 insertions(+), 59 deletions(-) create mode 100644 LiteDB.Benchmarks/Benchmarks/Spatial/SpatialQueryBenchmarks.cs create mode 100644 LiteDB.Benchmarks/Models/Spatial/SpatialDocument.cs create mode 100644 LiteDB.Benchmarks/Models/Spatial/SpatialDocumentGenerator.cs create mode 100644 LiteDB/Client/Mapper/Linq/TypeResolver/SpatialResolver.cs create mode 100644 LiteDB/Document/Expression/Methods/Spatial.cs create mode 100644 LiteDB/Spatial/SpatialMetadataStore.cs create mode 100644 LiteDB/Spatial/SpatialQueryBuilder.cs create mode 100644 docs/spatial-guide.md create mode 100644 samples/SpatialApiSample/Program.cs create mode 100644 samples/SpatialApiSample/SpatialApiSample.csproj diff --git a/LiteDB.Benchmarks/Benchmarks/Spatial/SpatialQueryBenchmarks.cs b/LiteDB.Benchmarks/Benchmarks/Spatial/SpatialQueryBenchmarks.cs new file mode 100644 index 000000000..a9633268a --- /dev/null +++ b/LiteDB.Benchmarks/Benchmarks/Spatial/SpatialQueryBenchmarks.cs @@ -0,0 +1,74 @@ +using System.Collections.Generic; +using System.IO; +using System.Linq; +using BenchmarkDotNet.Attributes; +using LiteDB.Benchmarks.Models.Spatial; +using LiteDB.Spatial; + +namespace LiteDB.Benchmarks.Benchmarks.Spatial +{ + [BenchmarkCategory(Constants.Categories.QUERIES)] + public class SpatialQueryBenchmarks : BenchmarkBase + { + private ILiteCollection _collection = null!; + private GeoPoint _center = null!; + private GeoPolygon _searchArea = null!; + private double _radiusMeters; + + [GlobalSetup] + public void GlobalSetup() + { + File.Delete(DatabasePath); + + DatabaseInstance = new LiteDatabase(ConnectionString()); + _collection = DatabaseInstance.GetCollection("places"); + + Spatial.EnsurePointIndex(_collection, x => x.Location); + Spatial.EnsureShapeIndex(_collection, x => x.Region); + Spatial.EnsureShapeIndex(_collection, x => x.Route); + + var documents = SpatialDocumentGenerator.Generate(DatasetSize); + _collection.Insert(documents); + + DatabaseInstance.Checkpoint(); + + _center = new GeoPoint(0, 0); + _radiusMeters = 25_000; + _searchArea = SpatialDocumentGenerator.BuildSearchPolygon(0, 0, 0.1); + } + + [Benchmark(Baseline = true)] + public List NearQuery() + { + return Spatial.Near(_collection, x => x.Location, _center, _radiusMeters).ToList(); + } + + [Benchmark] + public List BoundingBoxQuery() + { + return Spatial.WithinBoundingBox(_collection, x => x.Location, -0.2, -0.2, 0.2, 0.2).ToList(); + } + + [Benchmark] + public List PolygonContainmentQuery() + { + return Spatial.Within(_collection, x => x.Region, _searchArea).ToList(); + } + + [Benchmark] + public List RouteIntersectionQuery() + { + return Spatial.Intersects(_collection, x => x.Route, _searchArea).ToList(); + } + + [GlobalCleanup] + public void GlobalCleanup() + { + DatabaseInstance?.Checkpoint(); + DatabaseInstance?.Dispose(); + DatabaseInstance = null; + + File.Delete(DatabasePath); + } + } +} diff --git a/LiteDB.Benchmarks/Models/Spatial/SpatialDocument.cs b/LiteDB.Benchmarks/Models/Spatial/SpatialDocument.cs new file mode 100644 index 000000000..a0de19538 --- /dev/null +++ b/LiteDB.Benchmarks/Models/Spatial/SpatialDocument.cs @@ -0,0 +1,32 @@ +using LiteDB.Spatial; + +namespace LiteDB.Benchmarks.Models.Spatial +{ + public class SpatialDocument + { + public int Id { get; set; } + + public string Name { get; set; } = string.Empty; + + public GeoPoint Location { get; set; } = new GeoPoint(0, 0); + + public GeoPolygon Region { get; set; } = new GeoPolygon(new[] + { + new GeoPoint(0, 0), + new GeoPoint(0, 0.001), + new GeoPoint(0.001, 0.001), + new GeoPoint(0.001, 0), + new GeoPoint(0, 0) + }); + + public GeoLineString Route { get; set; } = new GeoLineString(new[] + { + new GeoPoint(0, 0), + new GeoPoint(0.001, 0.001) + }); + + internal long _gh { get; set; } + + internal double[] _mbb { get; set; } = Array.Empty(); + } +} diff --git a/LiteDB.Benchmarks/Models/Spatial/SpatialDocumentGenerator.cs b/LiteDB.Benchmarks/Models/Spatial/SpatialDocumentGenerator.cs new file mode 100644 index 000000000..6fc03f239 --- /dev/null +++ b/LiteDB.Benchmarks/Models/Spatial/SpatialDocumentGenerator.cs @@ -0,0 +1,78 @@ +using System; +using System.Collections.Generic; +using LiteDB.Spatial; + +namespace LiteDB.Benchmarks.Models.Spatial +{ + internal static class SpatialDocumentGenerator + { + public static List Generate(int count) + { + var random = new Random(1337); + var documents = new List(count); + + for (var i = 0; i < count; i++) + { + var lat = random.NextDouble() * 0.8 - 0.4; + var lon = random.NextDouble() * 0.8 - 0.4; + var location = new GeoPoint(lat, lon); + + var region = BuildSquare(location, random.NextDouble() * 0.05 + 0.01); + var route = BuildRoute(location, random); + + documents.Add(new SpatialDocument + { + Id = i + 1, + Name = $"Place #{i + 1}", + Location = location, + Region = region, + Route = route + }); + } + + return documents; + } + + public static GeoPolygon BuildSearchPolygon(double centerLat, double centerLon, double radiusDegrees) + { + var center = new GeoPoint(centerLat, centerLon); + return BuildSquare(center, radiusDegrees); + } + + private static GeoPolygon BuildSquare(GeoPoint center, double halfExtent) + { + var minLat = GeoMath.ClampLatitude(center.Lat - halfExtent); + var maxLat = GeoMath.ClampLatitude(center.Lat + halfExtent); + var minLon = GeoMath.NormalizeLongitude(center.Lon - halfExtent); + var maxLon = GeoMath.NormalizeLongitude(center.Lon + halfExtent); + + var points = new List + { + new GeoPoint(maxLat, minLon), + new GeoPoint(maxLat, maxLon), + new GeoPoint(minLat, maxLon), + new GeoPoint(minLat, minLon), + new GeoPoint(maxLat, minLon) + }; + + return new GeoPolygon(points); + } + + private static GeoLineString BuildRoute(GeoPoint start, Random random) + { + var midLat = start.Lat + random.NextDouble() * 0.1 - 0.05; + var midLon = start.Lon + random.NextDouble() * 0.1 - 0.05; + var endLat = start.Lat + random.NextDouble() * 0.2 - 0.1; + var endLon = start.Lon + random.NextDouble() * 0.2 - 0.1; + + var points = new List + { + start, + new GeoPoint(midLat, midLon), + new GeoPoint(endLat, endLon) + }; + + return new GeoLineString(points); + } + } +} diff --git a/LiteDB.Tests/Spatial/SpatialTests.cs b/LiteDB.Tests/Spatial/SpatialTests.cs index 797b2d170..b5b70362b 100644 --- a/LiteDB.Tests/Spatial/SpatialTests.cs +++ b/LiteDB.Tests/Spatial/SpatialTests.cs @@ -166,6 +166,73 @@ public void Intersects_Line_Touches_Polygon_Edge() Assert.NotEmpty(hits); } + [Fact] + public void EnsurePointIndex_Persists_Precision_Metadata() + { + using var db = new LiteDatabase(":memory:"); + var col = db.GetCollection("p"); + + SpatialApi.EnsurePointIndex(col, x => x.Location, 40); + + var meta = db.GetCollection("_spatial_meta").FindAll().ToList(); + + Assert.Single(meta); + Assert.Equal(40, meta[0]["precisionBits"].AsInt32); + } + + [Fact] + public void Linq_Near_Uses_Spatial_Operator() + { + using var db = new LiteDatabase(":memory:"); + var col = db.GetCollection("p"); + SpatialApi.EnsurePointIndex(col, x => x.Location); + + col.Insert(new[] + { + new Place { Name = "Center", Location = new GeoPoint(48.2082, 16.3738) }, + new Place { Name = "Far", Location = new GeoPoint(48.35, 16.7) } + }); + + var center = new GeoPoint(48.2082, 16.3738); + + var results = col.Query() + .Where(p => SpatialApi.Near(p.Location, center, 1_000)) + .ToList(); + + Assert.Single(results); + Assert.Equal("Center", results[0].Name); + } + + [Fact] + public void Linq_Within_Uses_Spatial_Operator() + { + using var db = new LiteDatabase(":memory:"); + var col = db.GetCollection("p"); + SpatialApi.EnsureShapeIndex(col, x => x.Location); + + var polygon = new GeoPolygon(new[] + { + new GeoPoint(48.25, 16.30), + new GeoPoint(48.25, 16.45), + new GeoPoint(48.15, 16.45), + new GeoPoint(48.15, 16.30), + new GeoPoint(48.25, 16.30) + }); + + col.Insert(new[] + { + new Place { Name = "Inside", Location = new GeoPoint(48.21, 16.37) }, + new Place { Name = "Outside", Location = new GeoPoint(48.30, 16.60) } + }); + + var results = col.Query() + .Where(p => SpatialApi.Within(p.Location, polygon)) + .ToList(); + + Assert.Single(results); + Assert.Equal("Inside", results[0].Name); + } + private class Place { public ObjectId Id { get; set; } diff --git a/LiteDB/Client/Mapper/Linq/LinqExpressionVisitor.cs b/LiteDB/Client/Mapper/Linq/LinqExpressionVisitor.cs index fe006012f..87fe84c38 100644 --- a/LiteDB/Client/Mapper/Linq/LinqExpressionVisitor.cs +++ b/LiteDB/Client/Mapper/Linq/LinqExpressionVisitor.cs @@ -30,7 +30,8 @@ internal class LinqExpressionVisitor : ExpressionVisitor [typeof(Regex)] = new RegexResolver(), [typeof(ObjectId)] = new ObjectIdResolver(), [typeof(String)] = new StringResolver(), - [typeof(Nullable)] = new NullableResolver() + [typeof(Nullable)] = new NullableResolver(), + [typeof(LiteDB.Spatial.Spatial)] = new SpatialResolver() }; private readonly BsonMapper _mapper; diff --git a/LiteDB/Client/Mapper/Linq/TypeResolver/SpatialResolver.cs b/LiteDB/Client/Mapper/Linq/TypeResolver/SpatialResolver.cs new file mode 100644 index 000000000..c25d75fef --- /dev/null +++ b/LiteDB/Client/Mapper/Linq/TypeResolver/SpatialResolver.cs @@ -0,0 +1,50 @@ +using System.Reflection; +using LiteDB.Spatial; +using SpatialMethods = LiteDB.Spatial.Spatial; + +namespace LiteDB +{ + internal class SpatialResolver : ITypeResolver + { + public string ResolveMethod(MethodInfo method) + { + if (method == null) + { + return null; + } + + if (method.DeclaringType != typeof(SpatialMethods)) + { + return null; + } + + var parameters = method.GetParameters(); + + if (method.Name == nameof(SpatialMethods.Near) && parameters.Length == 3 && parameters[0].ParameterType == typeof(GeoPoint)) + { + return "SPATIAL_NEAR(@0, @1, @2)"; + } + + if (method.Name == nameof(SpatialMethods.Within) && parameters.Length == 2 && parameters[0].ParameterType == typeof(GeoShape)) + { + return "SPATIAL_WITHIN(@0, @1)"; + } + + if (method.Name == nameof(SpatialMethods.Intersects) && parameters.Length == 2 && parameters[0].ParameterType == typeof(GeoShape)) + { + return "SPATIAL_INTERSECTS(@0, @1)"; + } + + if (method.Name == nameof(SpatialMethods.Contains) && parameters.Length == 2 && parameters[0].ParameterType == typeof(GeoShape)) + { + return "SPATIAL_CONTAINS_POINT(@0, @1)"; + } + + return null; + } + + public string ResolveMember(MemberInfo member) => null; + + public string ResolveCtor(ConstructorInfo ctor) => null; + } +} diff --git a/LiteDB/Document/Expression/Methods/Spatial.cs b/LiteDB/Document/Expression/Methods/Spatial.cs new file mode 100644 index 000000000..3be5bd916 --- /dev/null +++ b/LiteDB/Document/Expression/Methods/Spatial.cs @@ -0,0 +1,98 @@ +using LiteDB.Spatial; + +namespace LiteDB +{ + internal partial class BsonExpressionMethods + { + public static BsonValue SPATIAL_INTERSECTS_MBB(BsonValue mbb, BsonValue minLat, BsonValue minLon, BsonValue maxLat, BsonValue maxLon) + { + if (mbb == null || mbb.IsNull || !mbb.IsArray || mbb.AsArray.Count < 4) + { + return false; + } + + if (!minLat.IsNumber || !minLon.IsNumber || !maxLat.IsNumber || !maxLon.IsNumber) + { + return false; + } + + var candidate = new GeoBoundingBox(mbb.AsArray[0].AsDouble, mbb.AsArray[1].AsDouble, mbb.AsArray[2].AsDouble, mbb.AsArray[3].AsDouble); + var query = new GeoBoundingBox(minLat.AsDouble, minLon.AsDouble, maxLat.AsDouble, maxLon.AsDouble); + + return candidate.Intersects(query); + } + + public static BsonValue SPATIAL_NEAR(BsonValue candidate, BsonValue center, BsonValue radiusMeters) + { + if (!radiusMeters.IsNumber) + { + return false; + } + + var candidatePoint = ToGeoPoint(candidate); + var centerPoint = ToGeoPoint(center); + + if (candidatePoint == null || centerPoint == null) + { + return false; + } + + return Spatial.Spatial.Near(candidatePoint, centerPoint, radiusMeters.AsDouble); + } + + public static BsonValue SPATIAL_WITHIN(BsonValue candidate, BsonValue polygon) + { + var shape = ToGeoShape(candidate); + var area = ToGeoShape(polygon) as GeoPolygon; + + if (shape == null || area == null) + { + return false; + } + + return Spatial.Spatial.Within(shape, area); + } + + public static BsonValue SPATIAL_INTERSECTS(BsonValue candidate, BsonValue other) + { + var left = ToGeoShape(candidate); + var right = ToGeoShape(other); + + if (left == null || right == null) + { + return false; + } + + return Spatial.Spatial.Intersects(left, right); + } + + public static BsonValue SPATIAL_CONTAINS_POINT(BsonValue candidate, BsonValue point) + { + var shape = ToGeoShape(candidate); + var geoPoint = ToGeoPoint(point); + + if (shape == null || geoPoint == null) + { + return false; + } + + return Spatial.Spatial.Contains(shape, geoPoint); + } + + private static GeoShape ToGeoShape(BsonValue value) + { + if (value == null || value.IsNull || !value.IsDocument) + { + return null; + } + + return GeoJson.FromBson(value.AsDocument); + } + + private static GeoPoint ToGeoPoint(BsonValue value) + { + var shape = ToGeoShape(value); + return shape as GeoPoint; + } + } +} diff --git a/LiteDB/Spatial/GeoMath.cs b/LiteDB/Spatial/GeoMath.cs index 1554df9f0..2632f21ce 100644 --- a/LiteDB/Spatial/GeoMath.cs +++ b/LiteDB/Spatial/GeoMath.cs @@ -8,7 +8,7 @@ internal static class GeoMath private const double DegToRad = Math.PI / 180d; - internal const double EpsilonDegrees = 1e-9; + internal static double EpsilonDegrees => Spatial.Options.NumericToleranceDegrees; public static double ClampLatitude(double latitude) { diff --git a/LiteDB/Spatial/Spatial.cs b/LiteDB/Spatial/Spatial.cs index d6d103ba9..5411e715d 100644 --- a/LiteDB/Spatial/Spatial.cs +++ b/LiteDB/Spatial/Spatial.cs @@ -4,6 +4,7 @@ using System.Linq.Expressions; using System.Reflection; using LiteDB; +using LiteDB.Engine; namespace LiteDB.Spatial { @@ -14,7 +15,7 @@ public static class Spatial public static SpatialOptions Options { get; set; } = new SpatialOptions(); - public static void EnsurePointIndex(ILiteCollection collection, Expression> selector, int precisionBits = 52) + public static void EnsurePointIndex(ILiteCollection collection, Expression> selector, int precisionBits = 0) { if (collection == null) throw new ArgumentNullException(nameof(collection)); if (selector == null) throw new ArgumentNullException(nameof(selector)); @@ -23,6 +24,11 @@ public static void EnsurePointIndex(ILiteCollection collection, Expression var mapper = GetMapper(lite); EnsureMapperRegistration(mapper); + if (precisionBits <= 0) + { + precisionBits = Options.IndexPrecisionBits; + } + var getter = selector.Compile(); SpatialMapping.EnsureComputedMember(lite, "_gh", typeof(long), entity => @@ -40,6 +46,8 @@ public static void EnsurePointIndex(ILiteCollection collection, Expression SpatialMapping.EnsureBoundingBox(lite, entity => getter(entity)?.Normalize()); lite.EnsureIndex("_gh", BsonExpression.Create("$._gh")); + + SpatialMetadataStore.PersistPointIndexMetadata(lite, precisionBits); } public static void EnsureShapeIndex(ILiteCollection collection, Expression> selector) @@ -69,6 +77,73 @@ private static void EnsureShapeIndexInternal(ILiteCollection collection, F SpatialMapping.EnsureBoundingBox(lite, getter); } + public static bool Near(GeoPoint candidate, GeoPoint center, double radiusMeters) + { + if (center == null) throw new ArgumentNullException(nameof(center)); + if (radiusMeters < 0d) throw new ArgumentOutOfRangeException(nameof(radiusMeters)); + + if (candidate == null) + { + return false; + } + + var distance = GeoMath.DistanceMeters(center, candidate, Options.Distance); + return distance <= radiusMeters; + } + + public static bool Within(GeoShape candidate, GeoPolygon area) + { + if (area == null) throw new ArgumentNullException(nameof(area)); + + return candidate switch + { + GeoPoint point => Geometry.ContainsPoint(area, point), + GeoPolygon polygon => Geometry.Intersects(area, polygon) && polygon.Outer.All(p => Geometry.ContainsPoint(area, p)), + GeoLineString line => line.Points.All(p => Geometry.ContainsPoint(area, p)), + _ => false + }; + } + + public static bool Intersects(GeoShape candidate, GeoShape query) + { + if (candidate == null || query == null) + { + return false; + } + + return candidate switch + { + GeoPoint point when query is GeoPolygon polygon => Geometry.ContainsPoint(polygon, point), + GeoPoint point when query is GeoLineString line => Geometry.LineContainsPoint(line, point), + GeoPoint point when query is GeoPoint other => Math.Abs(point.Lat - other.Lat) < GeoMath.EpsilonDegrees && Math.Abs(point.Lon - other.Lon) < GeoMath.EpsilonDegrees, + GeoLineString line when query is GeoPolygon polygon => Geometry.Intersects(line, polygon), + GeoLineString line when query is GeoLineString other => Geometry.Intersects(line, other), + GeoLineString line when query is GeoPoint point => Geometry.LineContainsPoint(line, point), + GeoPolygon polygon when query is GeoPolygon other => Geometry.Intersects(polygon, other), + GeoPolygon polygon when query is GeoLineString line => Geometry.Intersects(line, polygon), + GeoPolygon polygon when query is GeoPoint point => Geometry.ContainsPoint(polygon, point), + _ => false + }; + } + + public static bool Contains(GeoShape candidate, GeoPoint point) + { + if (point == null) throw new ArgumentNullException(nameof(point)); + + if (candidate == null) + { + return false; + } + + return candidate switch + { + GeoPolygon polygon => Geometry.ContainsPoint(polygon, point), + GeoLineString line => Geometry.LineContainsPoint(line, point), + GeoPoint candidatePoint => Math.Abs(candidatePoint.Lat - point.Lat) < GeoMath.EpsilonDegrees && Math.Abs(candidatePoint.Lon - point.Lon) < GeoMath.EpsilonDegrees, + _ => false + }; + } + public static IEnumerable Near(ILiteCollection collection, Func selector, GeoPoint center, double radiusMeters, int? limit = null) { if (collection == null) throw new ArgumentNullException(nameof(collection)); @@ -80,12 +155,20 @@ public static IEnumerable Near(ILiteCollection collection, Func(); + var precisionBits = SpatialMetadataStore.GetPointIndexPrecision(lite); + var boundingBox = GeoMath.BoundingBoxForCircle(center, radiusMeters); + var ranges = SpatialIndexing.CoverBoundingBox(boundingBox, precisionBits, Options.MaxCoveringCells); + var rangePredicate = SpatialQueryBuilder.BuildRangePredicate(ranges); + var boundingPredicate = SpatialQueryBuilder.BuildBoundingBoxPredicate(boundingBox); + var predicate = SpatialQueryBuilder.CombineSpatialPredicates(rangePredicate, boundingPredicate); + + var source = predicate != null ? lite.Find(predicate) : lite.FindAll(); + var matches = new List<(T item, double distance)>(); - foreach (var item in lite.FindAll()) + foreach (var item in source) { var point = selector(item); - if (point == null) + if (point == null || !boundingBox.Contains(point)) { continue; } @@ -93,21 +176,23 @@ public static IEnumerable Near(ILiteCollection collection, Func x.distance.CompareTo(y.distance)); + matches.Sort((x, y) => x.distance.CompareTo(y.distance)); } + IEnumerable<(T item, double distance)> final = matches; + if (limit.HasValue) { - return candidates.Take(limit.Value).Select(x => x.item).ToList(); + final = matches.Take(limit.Value); } - return candidates.Select(x => x.item).ToList(); + return final.Select(x => x.item).ToList(); } public static IEnumerable WithinBoundingBox(ILiteCollection collection, Func selector, double minLat, double minLon, double maxLat, double maxLon) @@ -119,19 +204,20 @@ public static IEnumerable WithinBoundingBox(ILiteCollection collection, var mapper = GetMapper(lite); EnsureMapperRegistration(mapper); - var latitudeRange = (min: Math.Min(minLat, maxLat), max: Math.Max(minLat, maxLat)); - var longitudeRange = new LongitudeRange(minLon, maxLon); + var boundingBox = new GeoBoundingBox(minLat, minLon, maxLat, maxLon); + var precisionBits = SpatialMetadataStore.GetPointIndexPrecision(lite); + var ranges = SpatialIndexing.CoverBoundingBox(boundingBox, precisionBits, Options.MaxCoveringCells); + var rangePredicate = SpatialQueryBuilder.BuildRangePredicate(ranges); + var boundingPredicate = SpatialQueryBuilder.BuildBoundingBoxPredicate(boundingBox); + var predicate = SpatialQueryBuilder.CombineSpatialPredicates(rangePredicate, boundingPredicate); + + var source = predicate != null ? lite.Find(predicate) : lite.FindAll(); - return lite.FindAll() + return source .Where(entity => { var point = selector(entity); - if (point == null) - { - return false; - } - - return point.Lat >= latitudeRange.min && point.Lat <= latitudeRange.max && longitudeRange.Contains(point.Lon); + return point != null && boundingBox.Contains(point); }) .ToList(); } @@ -146,21 +232,14 @@ public static IEnumerable Within(ILiteCollection collection, Func + var boundingBox = area.GetBoundingBox(); + var boundingPredicate = SpatialQueryBuilder.BuildBoundingBoxPredicate(boundingBox); + var source = boundingPredicate != null ? lite.Find(boundingPredicate) : lite.FindAll(); + + return source.Where(entity => { var shape = selector(entity); - if (shape == null) - { - return false; - } - - return shape switch - { - GeoPoint point => Geometry.ContainsPoint(area, point), - GeoPolygon polygon => Geometry.Intersects(area, polygon) && polygon.Outer.All(p => Geometry.ContainsPoint(area, p)), - GeoLineString line => line.Points.All(p => Geometry.ContainsPoint(area, p)), - _ => false - }; + return shape != null && Within(shape, area); }).ToList(); } @@ -174,22 +253,14 @@ public static IEnumerable Intersects(ILiteCollection collection, Func + var boundingBox = query.GetBoundingBox(); + var boundingPredicate = SpatialQueryBuilder.BuildBoundingBoxPredicate(boundingBox); + var source = boundingPredicate != null ? lite.Find(boundingPredicate) : lite.FindAll(); + + return source.Where(entity => { var shape = selector(entity); - if (shape == null) - { - return false; - } - - return shape switch - { - GeoPoint point when query is GeoPolygon polygon => Geometry.ContainsPoint(polygon, point), - GeoLineString line when query is GeoPolygon polygon => Geometry.Intersects(line, polygon), - GeoPolygon polygon when query is GeoPolygon other => Geometry.Intersects(polygon, other), - GeoLineString line when query is GeoLineString other => Geometry.Intersects(line, other), - _ => false - }; + return shape != null && Intersects(shape, query); }).ToList(); } @@ -203,21 +274,14 @@ public static IEnumerable Contains(ILiteCollection collection, Func + var boundingBox = new GeoBoundingBox(point.Lat, point.Lon, point.Lat, point.Lon); + var boundingPredicate = SpatialQueryBuilder.BuildBoundingBoxPredicate(boundingBox); + var source = boundingPredicate != null ? lite.Find(boundingPredicate) : lite.FindAll(); + + return source.Where(entity => { var shape = selector(entity); - if (shape == null) - { - return false; - } - - return shape switch - { - GeoPolygon polygon => Geometry.ContainsPoint(polygon, point), - GeoLineString line => Geometry.LineContainsPoint(line, point), - GeoPoint candidate => Math.Abs(candidate.Lat - point.Lat) < GeoMath.EpsilonDegrees && Math.Abs(candidate.Lon - point.Lon) < GeoMath.EpsilonDegrees, - _ => false - }; + return shape != null && Contains(shape, point); }).ToList(); } @@ -231,6 +295,17 @@ private static LiteCollection GetLiteCollection(ILiteCollection collect throw new NotSupportedException("Spatial helpers require LiteCollection instances."); } + internal static ILiteEngine GetEngine(LiteCollection liteCollection) + { + var engineField = typeof(LiteCollection).GetField("_engine", BindingFlags.NonPublic | BindingFlags.Instance); + if (engineField == null) + { + return null; + } + + return (ILiteEngine)engineField.GetValue(liteCollection); + } + private static BsonMapper GetMapper(LiteCollection liteCollection) { var mapperField = typeof(LiteCollection).GetField("_mapper", BindingFlags.NonPublic | BindingFlags.Instance); diff --git a/LiteDB/Spatial/SpatialIndexing.cs b/LiteDB/Spatial/SpatialIndexing.cs index a021ba962..052e81220 100644 --- a/LiteDB/Spatial/SpatialIndexing.cs +++ b/LiteDB/Spatial/SpatialIndexing.cs @@ -1,4 +1,5 @@ using System; +using System.Collections.Generic; namespace LiteDB.Spatial { @@ -48,5 +49,106 @@ private static ulong Interleave(ulong x, ulong y, int bits) return result; } + + public static IReadOnlyList<(long Start, long End)> CoverBoundingBox(GeoBoundingBox box, int precisionBits, int maxCells) + { + if (precisionBits <= 0) + { + return Array.Empty<(long Start, long End)>(); + } + + maxCells = Math.Max(1, maxCells); + + var segments = SplitBoundingBox(box); + var ranges = new List<(long Start, long End)>(); + + var latCells = Math.Max(1, (int)Math.Round(Math.Sqrt(maxCells))); + var lonCells = Math.Max(1, maxCells / latCells); + + foreach (var segment in segments) + { + var latStep = (segment.MaxLat - segment.MinLat) / latCells; + var lonStep = (segment.MaxLon - segment.MinLon) / lonCells; + + if (latStep == 0) + { + latStep = segment.MaxLat - segment.MinLat; + } + + if (lonStep == 0) + { + lonStep = segment.MaxLon - segment.MinLon; + } + + for (var latIndex = 0; latIndex < latCells; latIndex++) + { + var minLat = segment.MinLat + latStep * latIndex; + var maxLat = latIndex == latCells - 1 ? segment.MaxLat : minLat + latStep; + + for (var lonIndex = 0; lonIndex < lonCells; lonIndex++) + { + var minLon = segment.MinLon + lonStep * lonIndex; + var maxLon = lonIndex == lonCells - 1 ? segment.MaxLon : minLon + lonStep; + + var start = ComputeMorton(new GeoPoint(minLat, minLon), precisionBits); + var end = ComputeMorton(new GeoPoint(maxLat, maxLon), precisionBits); + + if (start > end) + { + (start, end) = (end, start); + } + + ranges.Add((start, end)); + } + } + } + + return MergeRanges(ranges); + } + + public static IReadOnlyList SplitBoundingBox(GeoBoundingBox box) + { + if (box.MaxLon >= box.MinLon) + { + return new[] { box }; + } + + var first = new GeoBoundingBox(box.MinLat, box.MinLon, box.MaxLat, 180d); + var second = new GeoBoundingBox(box.MinLat, -180d, box.MaxLat, box.MaxLon); + + return new[] { first, second }; + } + + private static IReadOnlyList<(long Start, long End)> MergeRanges(List<(long Start, long End)> ranges) + { + if (ranges.Count == 0) + { + return ranges; + } + + ranges.Sort((a, b) => a.Start.CompareTo(b.Start)); + + var merged = new List<(long Start, long End)>(ranges.Count); + var current = ranges[0]; + + for (var i = 1; i < ranges.Count; i++) + { + var next = ranges[i]; + + if (next.Start <= current.End + 1) + { + current = (current.Start, Math.Max(current.End, next.End)); + } + else + { + merged.Add(current); + current = next; + } + } + + merged.Add(current); + + return merged; + } } } diff --git a/LiteDB/Spatial/SpatialMetadataStore.cs b/LiteDB/Spatial/SpatialMetadataStore.cs new file mode 100644 index 000000000..977adfe9c --- /dev/null +++ b/LiteDB/Spatial/SpatialMetadataStore.cs @@ -0,0 +1,129 @@ +using System; +using System.Collections.Concurrent; +using System.Runtime.CompilerServices; +using LiteDB.Engine; +using LiteDB; + +namespace LiteDB.Spatial +{ + internal static class SpatialMetadataStore + { + private const string MetadataCollection = "_spatial_meta"; + + private static readonly ConcurrentDictionary _cache = new(); + + public static void PersistPointIndexMetadata(LiteCollection collection, int precisionBits) + { + var engine = Spatial.GetEngine(collection); + if (engine == null) + { + return; + } + + var key = EngineCollectionKey.Create(engine, collection.Name); + + var document = new BsonDocument + { + ["_id"] = new BsonValue($"{collection.Name}._gh"), + ["collection"] = collection.Name, + ["index"] = "_gh", + ["precisionBits"] = precisionBits, + ["updatedUtc"] = DateTime.UtcNow + }; + + engine.Upsert(MetadataCollection, new[] { document }, BsonAutoId.ObjectId); + + _cache[key] = new SpatialIndexMetadata(precisionBits); + } + + public static int GetPointIndexPrecision(LiteCollection collection) + { + var engine = Spatial.GetEngine(collection); + if (engine == null) + { + return Spatial.Options.IndexPrecisionBits; + } + + var key = EngineCollectionKey.Create(engine, collection.Name); + + if (_cache.TryGetValue(key, out var metadata)) + { + return metadata.PrecisionBits; + } + + var predicate = Query.And( + Query.EQ("collection", new BsonValue(collection.Name)), + Query.EQ("index", new BsonValue("_gh"))); + + try + { + using var reader = engine.Query(MetadataCollection, new Query { Where = { predicate } }); + + if (reader.Read()) + { + var document = reader.Current.AsDocument; + + if (document.TryGetValue("precisionBits", out var precisionValue) && precisionValue.IsInt32) + { + var precisionBits = precisionValue.AsInt32; + _cache[key] = new SpatialIndexMetadata(precisionBits); + return precisionBits; + } + } + } + catch (LiteException) + { + // Metadata collection may not exist yet; fall back to options. + } + + return Spatial.Options.IndexPrecisionBits; + } + + private readonly struct SpatialIndexMetadata + { + public SpatialIndexMetadata(int precisionBits) + { + this.PrecisionBits = precisionBits; + } + + public int PrecisionBits { get; } + } + + private readonly struct EngineCollectionKey : IEquatable + { + private EngineCollectionKey(ILiteEngine engine, string collection) + { + _engine = engine; + _collection = collection; + } + + private readonly ILiteEngine _engine; + private readonly string _collection; + + public static EngineCollectionKey Create(ILiteEngine engine, string collection) + { + return new EngineCollectionKey(engine, collection); + } + + public bool Equals(EngineCollectionKey other) + { + return ReferenceEquals(_engine, other._engine) && string.Equals(_collection, other._collection, StringComparison.OrdinalIgnoreCase); + } + + public override bool Equals(object obj) + { + return obj is EngineCollectionKey other && Equals(other); + } + + public override int GetHashCode() + { + unchecked + { + var hash = RuntimeHelpers.GetHashCode(_engine); + hash = (hash * 397) ^ (_collection == null ? 0 : StringComparer.OrdinalIgnoreCase.GetHashCode(_collection)); + return hash; + } + } + } + } +} diff --git a/LiteDB/Spatial/SpatialOptions.cs b/LiteDB/Spatial/SpatialOptions.cs index 8f5e3278d..802871d75 100644 --- a/LiteDB/Spatial/SpatialOptions.cs +++ b/LiteDB/Spatial/SpatialOptions.cs @@ -21,5 +21,9 @@ public sealed class SpatialOptions public int MaxCoveringCells { get; set; } = 32; public AngleUnit AngleUnit { get; set; } = AngleUnit.Degrees; + + public int IndexPrecisionBits { get; set; } = 52; + + public double NumericToleranceDegrees { get; set; } = 1e-9; } } diff --git a/LiteDB/Spatial/SpatialQueryBuilder.cs b/LiteDB/Spatial/SpatialQueryBuilder.cs new file mode 100644 index 000000000..1a1916a69 --- /dev/null +++ b/LiteDB/Spatial/SpatialQueryBuilder.cs @@ -0,0 +1,94 @@ +using System.Collections.Generic; +using LiteDB; + +namespace LiteDB.Spatial +{ + internal static class SpatialQueryBuilder + { + public static BsonExpression BuildRangePredicate(IReadOnlyList<(long Start, long End)> ranges) + { + if (ranges == null || ranges.Count == 0) + { + return null; + } + + var expressions = new List(ranges.Count); + + foreach (var range in ranges) + { + expressions.Add(Query.Between("$._gh", new BsonValue(range.Start), new BsonValue(range.End))); + } + + return CombineOr(expressions); + } + + public static BsonExpression BuildBoundingBoxPredicate(GeoBoundingBox box) + { + if (box.MaxLon < box.MinLon) + { + return null; + } + + var parameters = new[] + { + new BsonValue(box.MaxLat), + new BsonValue(box.MinLat), + new BsonValue(box.MaxLon), + new BsonValue(box.MinLon) + }; + + const string predicate = "($._mbb != null) AND $._mbb[0] <= @0 AND $._mbb[2] >= @1 AND $._mbb[1] <= @2 AND $._mbb[3] >= @3"; + + return BsonExpression.Create(predicate, parameters); + } + + public static BsonExpression CombineSpatialPredicates(BsonExpression rangeExpression, BsonExpression boundingExpression) + { + if (rangeExpression == null) + { + return boundingExpression; + } + + if (boundingExpression == null) + { + return rangeExpression; + } + + var parameters = new BsonDocument(); + + if (boundingExpression.Parameters != null) + { + foreach (var parameter in boundingExpression.Parameters) + { + parameters[parameter.Key] = parameter.Value; + } + } + + var source = $"({rangeExpression.Source}) AND ({boundingExpression.Source})"; + + return BsonExpression.Create(source, parameters); + } + + private static BsonExpression CombineOr(IReadOnlyList expressions) + { + if (expressions.Count == 0) + { + return null; + } + + if (expressions.Count == 1) + { + return expressions[0]; + } + + var buffer = new BsonExpression[expressions.Count]; + + for (var i = 0; i < expressions.Count; i++) + { + buffer[i] = expressions[i]; + } + + return Query.Or(buffer); + } + } +} diff --git a/docs/spatial-guide.md b/docs/spatial-guide.md new file mode 100644 index 000000000..4cb0420af --- /dev/null +++ b/docs/spatial-guide.md @@ -0,0 +1,129 @@ +# Spatial Indexing Guide + +LiteDB's spatial tooling now supports index-aware queries, expression operators, and ready-to-run samples. This guide walks through enabling spatial indexes, querying data, and composing applications that take advantage of the new API surface. + +## 1. Preparing Collections + +```csharp +using LiteDB; +using LiteDB.Spatial; + +using var db = new LiteDatabase("Filename=geo.db;Mode=Shared"); +var places = db.GetCollection("places"); + +// Persist precision metadata and computed members +Spatial.EnsurePointIndex(places, x => x.Location); +Spatial.EnsureShapeIndex(places, x => x.Footprint); +``` + +`EnsurePointIndex` now records the Morton precision that was used so future queries can translate shapes into the correct `_gh` windows without additional configuration. + +```csharp +public class Place +{ + public int Id { get; set; } + public string Name { get; set; } = string.Empty; + public GeoPoint Location { get; set; } = new GeoPoint(0, 0); + public GeoPolygon Footprint { get; set; } = SquareAround(0, 0, 0.05); + + internal long _gh { get; set; } + internal double[] _mbb { get; set; } = Array.Empty(); +} + +static GeoPolygon SquareAround(double lat, double lon, double halfExtent) +{ + var topLeft = new GeoPoint(lat + halfExtent, lon - halfExtent); + var topRight = new GeoPoint(lat + halfExtent, lon + halfExtent); + var bottomRight = new GeoPoint(lat - halfExtent, lon + halfExtent); + var bottomLeft = new GeoPoint(lat - halfExtent, lon - halfExtent); + + return new GeoPolygon(new[] { topLeft, topRight, bottomRight, bottomLeft, topLeft }); +} +``` + +## 2. Index-Aware Queries + +### Radius Searches + +`Spatial.Near` now projects circle queries into Morton range scans and `_mbb` filters before falling back to geometry checks. This avoids `FindAll()` enumeration even on large collections. + +```csharp +var vienna = new GeoPoint(48.2082, 16.3738); +var withinFiveKm = Spatial.Near(places, x => x.Location, vienna, radiusMeters: 5_000).ToList(); +``` + +### Bounding Boxes + +Bounding-box queries reuse the same range generator and anti-meridian aware filters: + +```csharp +var hits = Spatial.WithinBoundingBox(places, x => x.Location, 47.9, 16.1, 48.4, 16.6).ToList(); +``` + +### Polygon Containment & Intersections + +Shape-based queries can rely on the lightweight `_mbb` predicate that gets folded into the pipeline before precise geometry calculations: + +```csharp +var downtown = SquareAround(48.2082, 16.3738, 0.15); +var inside = Spatial.Within(places, x => x.Footprint, downtown).ToList(); +``` + +## 3. Expression & LINQ Operators + +Spatial functions participate in the LINQ translator and the expression engine via the new operators: + +| C# Call | Bson Expression | +| --- | --- | +| `Spatial.Near(doc.Location, center, 1000)` | `SPATIAL_NEAR($.Location, @0, @1)` | +| `Spatial.Within(doc.Footprint, polygon)` | `SPATIAL_WITHIN($.Footprint, @0)` | +| `Spatial.Intersects(doc.Route, query)` | `SPATIAL_INTERSECTS($.Route, @0)` | +| `Spatial.Contains(doc.Footprint, point)` | `SPATIAL_CONTAINS_POINT($.Footprint, @0)` | + +```csharp +var linq = places.Query() + .Where(p => Spatial.Near(p.Location, vienna, 2_000)) + .Select(p => new { p.Name, p.Location }) + .ToList(); +``` + +Each operator pairs with the indexed `_gh` ranges captured by `EnsurePointIndex`, maintaining fast candidate pruning. + +## 4. Benchmarking Spatial Pipelines + +`LiteDB.Benchmarks` ships with a `SpatialQueryBenchmarks` suite that tracks radius, bounding-box, containment, and intersection workloads across dataset sizes. Run: + +```bash +dotnet run --project LiteDB.Benchmarks -c Release --filter "SpatialQueryBenchmarks" +``` + +The new benchmarks emit allocations and wall-clock metrics so regressions are easy to spot as spatial features evolve. + +## 5. Sample REST API + +A minimal API showcasing radius and polygon queries lives under `samples/SpatialApiSample`. Seed and query via: + +```bash +dotnet run --project samples/SpatialApiSample +# In another terminal +curl -X POST http://localhost:5000/seed +curl "http://localhost:5000/places/near?lat=48.2&lon=16.37&radiusKm=5" +``` + +The endpoint reuses the shared helpers, ensuring metadata is created automatically and results arrive sorted by distance. + +## 6. Troubleshooting & Options + +`SpatialOptions` exposes tunables for query precision and numeric tolerances: + +```csharp +Spatial.Options = new SpatialOptions +{ + IndexPrecisionBits = 48, + NumericToleranceDegrees = 1e-8, + MaxCoveringCells = 64, + Distance = DistanceFormula.Vincenty +}; +``` + +Changing `IndexPrecisionBits` updates persisted metadata the next time `EnsurePointIndex` runs, so the engine always knows how to slice query ranges. Adjust `NumericToleranceDegrees` if your datasets require more relaxed comparisons for noisy coordinates. diff --git a/docs/spatial-roadmap.md b/docs/spatial-roadmap.md index 8266c80ad..8256a2e04 100644 --- a/docs/spatial-roadmap.md +++ b/docs/spatial-roadmap.md @@ -39,10 +39,12 @@ dotnet test LiteDB.Tests/LiteDB.Tests.csproj -f net8.0 --filter FullyQualifiedNa - Translate query shapes (circles, polygons) into `_gh` range windows and `_mbb` filters. - Hook range scans into the query pipeline to avoid `FindAll()` enumeration. - Benchmark on large datasets to confirm IO/CPU gains. +- ✅ Implemented in `SpatialIndexing.CoverBoundingBox`, `SpatialQueryBuilder`, and the updated `Spatial.Near`/`Within*` helpers. Range windows respect persisted precision metadata and combine with `_mbb` predicates before geometry checks. ### 2. LINQ & Expression Support - Introduce spatial operators into the BsonExpression engine (e.g., `$near`, `$within`, `$intersects`). - Extend the LINQ translator to recognise Spatial methods and emit the new operators. +- ✅ Added `SPATIAL_NEAR`, `SPATIAL_WITHIN`, `SPATIAL_INTERSECTS`, and `SPATIAL_CONTAINS_POINT` methods alongside a `SpatialResolver` so `Spatial.*` calls translate directly inside `LiteCollection.Query()` pipelines. ### 3. 3D & Extended Geometry - Design GeoPoint3D / GeoBoundingBox3D / GeoPolyhedron structures. @@ -52,6 +54,7 @@ dotnet test LiteDB.Tests/LiteDB.Tests.csproj -f net8.0 --filter FullyQualifiedNa ### 4. Precision & Options - Surface defaults via SpatialOptions (index precision, tolerance, distance formula). - Persist precision metadata alongside index definitions for smarter range calculations. +- ✅ `SpatialOptions` now exposes `IndexPrecisionBits` and `NumericToleranceDegrees`; `EnsurePointIndex` writes metadata into `_spatial_meta` and reuses it during query planning. ### 5. Migration & Tooling - Offer shell commands or utility APIs to backfill `_gh`/`_mbb` for existing datasets. @@ -60,10 +63,12 @@ dotnet test LiteDB.Tests/LiteDB.Tests.csproj -f net8.0 --filter FullyQualifiedNa ### 6. Performance Tracking - Add BenchmarkDotNet scenarios targeting Near/Within/Intersects across dataset sizes. - Monitor allocations and wall-clock time before/after index-aware implementation. +- ✅ `SpatialQueryBenchmarks` exercises radius, bounding-box, containment, and intersection workloads to capture allocations and timings for the new pipeline. ### 7. Documentation & Samples - Expand docs with tutorials covering spatial CRUD, indexing, and querying patterns. - Provide sample apps (e.g., REST endpoint performing radius searches). +- ✅ `docs/spatial-guide.md` consolidates tutorials, and `samples/SpatialApiSample` offers a Minimal API demonstrating seeding, radius searches, and polygon filters. ## Open Questions - Should `_gh`/`_mbb` remain internal fields or become part of the public query DSL? diff --git a/samples/SpatialApiSample/Program.cs b/samples/SpatialApiSample/Program.cs new file mode 100644 index 000000000..3bf760f54 --- /dev/null +++ b/samples/SpatialApiSample/Program.cs @@ -0,0 +1,103 @@ +using LiteDB; +using LiteDB.Spatial; + +var builder = WebApplication.CreateBuilder(args); +var app = builder.Build(); + +const string DatabasePath = "spatial-sample.db"; + +app.MapPost("/seed", () => +{ + using var db = new LiteDatabase(DatabasePath); + var places = db.GetCollection("places"); + Spatial.EnsurePointIndex(places, x => x.Location); + Spatial.EnsureShapeIndex(places, x => x.Coverage); + + if (places.Count() > 0) + { + return Results.Ok(new { message = "Database already seeded." }); + } + + var vienna = new Place + { + Name = "Vienna", + Location = new GeoPoint(48.2082, 16.3738), + Coverage = SquareAround(48.2082, 16.3738, 0.2) + }; + + var bratislava = new Place + { + Name = "Bratislava", + Location = new GeoPoint(48.1486, 17.1077), + Coverage = SquareAround(48.1486, 17.1077, 0.15) + }; + + places.Insert(new[] { vienna, bratislava }); + + return Results.Ok(new { message = "Seeded" }); +}); + +app.MapGet("/places/near", (double lat, double lon, double radiusKm) => +{ + using var db = new LiteDatabase(DatabasePath); + var places = db.GetCollection("places"); + Spatial.EnsurePointIndex(places, x => x.Location); + + var center = new GeoPoint(lat, lon); + var radiusMeters = radiusKm * 1000; + + var results = Spatial.Near(places, x => x.Location, center, radiusMeters) + .Select(x => new { x.Name, x.Location.Lat, x.Location.Lon }) + .ToList(); + + return Results.Ok(results); +}); + +app.MapGet("/places/within", () => +{ + using var db = new LiteDatabase(DatabasePath); + var places = db.GetCollection("places"); + Spatial.EnsureShapeIndex(places, x => x.Coverage); + + var polygon = SquareAround(48.2, 16.35, 0.25); + + var results = Spatial.Within(places, x => x.Coverage, polygon) + .Select(x => new { x.Name }) + .ToList(); + + return Results.Ok(results); +}); + +app.Run(); + +static GeoPolygon SquareAround(double lat, double lon, double halfExtent) +{ + var topLeft = new GeoPoint(lat + halfExtent, lon - halfExtent); + var topRight = new GeoPoint(lat + halfExtent, lon + halfExtent); + var bottomRight = new GeoPoint(lat - halfExtent, lon + halfExtent); + var bottomLeft = new GeoPoint(lat - halfExtent, lon - halfExtent); + + return new GeoPolygon(new[] + { + topLeft, + topRight, + bottomRight, + bottomLeft, + topLeft + }); +} + +public class Place +{ + public int Id { get; set; } + + public string Name { get; set; } = string.Empty; + + public GeoPoint Location { get; set; } = new GeoPoint(0, 0); + + public GeoPolygon Coverage { get; set; } = SquareAround(0, 0, 0.1); + + internal long _gh { get; set; } + + internal double[] _mbb { get; set; } = Array.Empty(); +} diff --git a/samples/SpatialApiSample/SpatialApiSample.csproj b/samples/SpatialApiSample/SpatialApiSample.csproj new file mode 100644 index 000000000..5818ac390 --- /dev/null +++ b/samples/SpatialApiSample/SpatialApiSample.csproj @@ -0,0 +1,10 @@ + + + net8.0 + enable + enable + + + + + From 6d185ad34385a063179e9686607a1a4078f59621 Mon Sep 17 00:00:00 2001 From: JKamsker <11245306+JKamsker@users.noreply.github.com> Date: Mon, 6 Oct 2025 01:09:54 +0200 Subject: [PATCH 4/8] Enhance queryoptimizations --- LiteDB/Engine/Query/QueryOptimization.cs | 28 ++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/LiteDB/Engine/Query/QueryOptimization.cs b/LiteDB/Engine/Query/QueryOptimization.cs index 92b013997..1cae2be3c 100644 --- a/LiteDB/Engine/Query/QueryOptimization.cs +++ b/LiteDB/Engine/Query/QueryOptimization.cs @@ -85,7 +85,7 @@ void add(BsonExpression predicate) } // add expression in where list breaking AND statments - if (predicate.IsPredicate || predicate.Type == BsonExpressionType.Or) + if (predicate.IsPredicate || predicate.Type == BsonExpressionType.Or || IsSpatialPredicate(predicate)) { _terms.Add(predicate); } @@ -110,6 +110,30 @@ void add(BsonExpression predicate) } } + private static bool IsSpatialPredicate(BsonExpression expression) + { + if (expression == null || expression.Type != BsonExpressionType.Call) + { + return false; + } + + var source = expression.Source?.Trim() ?? string.Empty; + var nameEnd = source.IndexOf('('); + var name = nameEnd >= 0 ? source.Substring(0, nameEnd) : source; + + switch (name.ToUpperInvariant()) + { + case "SPATIAL_INTERSECTS": + case "SPATIAL_INTERSECTS_MBB": + case "SPATIAL_CONTAINS_POINT": + case "SPATIAL_NEAR": + case "SPATIAL_WITHIN": + return true; + default: + return false; + } + } + /// /// Do some pre-defined optimization on terms to convert expensive filter in indexable filter /// @@ -383,4 +407,4 @@ private void DefineIncludes() } } } -} \ No newline at end of file +} From 29ccc6f73bfcebd8a520da165042807a4d4b9458 Mon Sep 17 00:00:00 2001 From: Jonas Kamsker <11245306+JKamsker@users.noreply.github.com> Date: Mon, 6 Oct 2025 07:45:57 +0200 Subject: [PATCH 5/8] Add spatial tolerance handling and expressions --- .../Linq/TypeResolver/SpatialResolver.cs | 58 ++++++++- LiteDB/Document/Expression/Methods/Spatial.cs | 121 ++++++++++++++---- LiteDB/Spatial/GeoBoundingBox.cs | 18 +++ LiteDB/Spatial/GeoMath.cs | 2 +- LiteDB/Spatial/Spatial.cs | 106 +++++++-------- LiteDB/Spatial/SpatialExpressions.cs | 95 ++++++++++++++ LiteDB/Spatial/SpatialMetadataStore.cs | 4 +- LiteDB/Spatial/SpatialOptions.cs | 31 ++++- 8 files changed, 350 insertions(+), 85 deletions(-) create mode 100644 LiteDB/Spatial/SpatialExpressions.cs diff --git a/LiteDB/Client/Mapper/Linq/TypeResolver/SpatialResolver.cs b/LiteDB/Client/Mapper/Linq/TypeResolver/SpatialResolver.cs index c25d75fef..fe8ac08e3 100644 --- a/LiteDB/Client/Mapper/Linq/TypeResolver/SpatialResolver.cs +++ b/LiteDB/Client/Mapper/Linq/TypeResolver/SpatialResolver.cs @@ -1,3 +1,4 @@ +using System; using System.Reflection; using LiteDB.Spatial; using SpatialMethods = LiteDB.Spatial.Spatial; @@ -13,16 +14,49 @@ public string ResolveMethod(MethodInfo method) return null; } - if (method.DeclaringType != typeof(SpatialMethods)) + if (method.DeclaringType == typeof(SpatialExpressions)) { - return null; + return ResolveSpatialExpressions(method); + } + + if (method.DeclaringType == typeof(SpatialMethods)) + { + return ResolveSpatialMethods(method); + } + + return null; + } + + public string ResolveMember(MemberInfo member) => null; + + public string ResolveCtor(ConstructorInfo ctor) => null; + + private static string ResolveSpatialExpressions(MethodInfo method) + { + switch (method.Name) + { + case nameof(SpatialExpressions.Near): + return ResolveNearPattern(method); + case nameof(SpatialExpressions.Within): + return "SPATIAL_WITHIN(@0, @1)"; + case nameof(SpatialExpressions.Intersects): + return "SPATIAL_INTERSECTS(@0, @1)"; + case nameof(SpatialExpressions.Contains): + return "SPATIAL_CONTAINS(@0, @1)"; + case nameof(SpatialExpressions.WithinBoundingBox): + return "SPATIAL_WITHIN_BOX(@0, @1, @2, @3, @4)"; } + return null; + } + + private static string ResolveSpatialMethods(MethodInfo method) + { var parameters = method.GetParameters(); if (method.Name == nameof(SpatialMethods.Near) && parameters.Length == 3 && parameters[0].ParameterType == typeof(GeoPoint)) { - return "SPATIAL_NEAR(@0, @1, @2)"; + return ResolveNearPattern(method); } if (method.Name == nameof(SpatialMethods.Within) && parameters.Length == 2 && parameters[0].ParameterType == typeof(GeoShape)) @@ -43,8 +77,22 @@ public string ResolveMethod(MethodInfo method) return null; } - public string ResolveMember(MemberInfo member) => null; + private static string ResolveNearPattern(MethodInfo method) + { + var parameters = method.GetParameters(); - public string ResolveCtor(ConstructorInfo ctor) => null; + if (parameters.Length == 3) + { + var formula = SpatialMethods.Options.Distance.ToString(); + return $"SPATIAL_NEAR(@0, @1, @2, '{formula}')"; + } + + if (parameters.Length == 4) + { + return "SPATIAL_NEAR(@0, @1, @2, @3)"; + } + + throw new NotSupportedException("Unsupported overload for spatial Near expression"); + } } } diff --git a/LiteDB/Document/Expression/Methods/Spatial.cs b/LiteDB/Document/Expression/Methods/Spatial.cs index 3be5bd916..1039c6af3 100644 --- a/LiteDB/Document/Expression/Methods/Spatial.cs +++ b/LiteDB/Document/Expression/Methods/Spatial.cs @@ -1,3 +1,4 @@ +using System; using LiteDB.Spatial; namespace LiteDB @@ -6,7 +7,7 @@ internal partial class BsonExpressionMethods { public static BsonValue SPATIAL_INTERSECTS_MBB(BsonValue mbb, BsonValue minLat, BsonValue minLon, BsonValue maxLat, BsonValue maxLon) { - if (mbb == null || mbb.IsNull || !mbb.IsArray || mbb.AsArray.Count < 4) + if (mbb == null || mbb.IsNull || !mbb.IsArray || mbb.AsArray.Count != 4) { return false; } @@ -16,83 +17,157 @@ public static BsonValue SPATIAL_INTERSECTS_MBB(BsonValue mbb, BsonValue minLat, return false; } - var candidate = new GeoBoundingBox(mbb.AsArray[0].AsDouble, mbb.AsArray[1].AsDouble, mbb.AsArray[2].AsDouble, mbb.AsArray[3].AsDouble); + var candidate = ToBoundingBox(mbb); var query = new GeoBoundingBox(minLat.AsDouble, minLon.AsDouble, maxLat.AsDouble, maxLon.AsDouble); return candidate.Intersects(query); } - public static BsonValue SPATIAL_NEAR(BsonValue candidate, BsonValue center, BsonValue radiusMeters) + public static BsonValue SPATIAL_MBB_INTERSECTS(BsonValue mbb, BsonValue minLat, BsonValue minLon, BsonValue maxLat, BsonValue maxLon) { - if (!radiusMeters.IsNumber) + return SPATIAL_INTERSECTS_MBB(mbb, minLat, minLon, maxLat, maxLon); + } + + public static BsonValue SPATIAL_NEAR(BsonValue candidate, BsonValue center, BsonValue radius) + { + return SPATIAL_NEAR(candidate, center, radius, BsonValue.Null); + } + + public static BsonValue SPATIAL_NEAR(BsonValue candidate, BsonValue center, BsonValue radius, BsonValue formula) + { + if (!radius.IsNumber) { return false; } - var candidatePoint = ToGeoPoint(candidate); - var centerPoint = ToGeoPoint(center); + var candidatePoint = ToShape(candidate) as GeoPoint; + var centerPoint = ToShape(center) as GeoPoint; if (candidatePoint == null || centerPoint == null) { return false; } - return Spatial.Spatial.Near(candidatePoint, centerPoint, radiusMeters.AsDouble); + var distanceFormula = ParseFormula(formula); + var radiusMeters = radius.AsDouble; + + return SpatialExpressions.Near(candidatePoint, centerPoint, radiusMeters, distanceFormula); + } + + public static BsonValue SPATIAL_WITHIN_BOX(BsonValue value, BsonValue minLat, BsonValue minLon, BsonValue maxLat, BsonValue maxLon) + { + var shape = ToShape(value); + if (shape == null) + { + return false; + } + + var box = new GeoBoundingBox(minLat.AsDouble, minLon.AsDouble, maxLat.AsDouble, maxLon.AsDouble); + + return shape switch + { + GeoPoint point => SpatialExpressions.WithinBoundingBox(point, box.MinLat, box.MinLon, box.MaxLat, box.MaxLon), + GeoShape geoShape => geoShape.GetBoundingBox().Intersects(box), + _ => false + }; } public static BsonValue SPATIAL_WITHIN(BsonValue candidate, BsonValue polygon) { - var shape = ToGeoShape(candidate); - var area = ToGeoShape(polygon) as GeoPolygon; + var shape = ToShape(candidate); + var area = ToShape(polygon) as GeoPolygon; if (shape == null || area == null) { return false; } - return Spatial.Spatial.Within(shape, area); + return SpatialExpressions.Within(shape, area); } public static BsonValue SPATIAL_INTERSECTS(BsonValue candidate, BsonValue other) { - var left = ToGeoShape(candidate); - var right = ToGeoShape(other); + var left = ToShape(candidate); + var right = ToShape(other); if (left == null || right == null) { return false; } - return Spatial.Spatial.Intersects(left, right); + return SpatialExpressions.Intersects(left, right); } - public static BsonValue SPATIAL_CONTAINS_POINT(BsonValue candidate, BsonValue point) + public static BsonValue SPATIAL_CONTAINS(BsonValue candidate, BsonValue point) { - var shape = ToGeoShape(candidate); - var geoPoint = ToGeoPoint(point); + var shape = ToShape(candidate); + var geoPoint = ToShape(point) as GeoPoint; if (shape == null || geoPoint == null) { return false; } - return Spatial.Spatial.Contains(shape, geoPoint); + return SpatialExpressions.Contains(shape, geoPoint); + } + + public static BsonValue SPATIAL_CONTAINS_POINT(BsonValue candidate, BsonValue point) + { + return SPATIAL_CONTAINS(candidate, point); + } + + private static GeoBoundingBox ToBoundingBox(BsonValue value) + { + var array = value.AsArray; + return new GeoBoundingBox(array[0].AsDouble, array[1].AsDouble, array[2].AsDouble, array[3].AsDouble); } - private static GeoShape ToGeoShape(BsonValue value) + private static GeoShape ToShape(BsonValue value) { - if (value == null || value.IsNull || !value.IsDocument) + if (value == null || value.IsNull) { return null; } - return GeoJson.FromBson(value.AsDocument); + if (value.IsDocument) + { + return GeoJson.FromBson(value.AsDocument); + } + + if (value.IsArray && value.AsArray.Count >= 2) + { + var array = value.AsArray; + var lon = array[0].AsDouble; + var lat = array[1].AsDouble; + return new GeoPoint(lat, lon); + } + + if (value.RawValue is GeoShape shape) + { + return shape; + } + + return null; } - private static GeoPoint ToGeoPoint(BsonValue value) + private static DistanceFormula ParseFormula(BsonValue formula) { - var shape = ToGeoShape(value); - return shape as GeoPoint; + if (formula == null || formula.IsNull) + { + return Spatial.Spatial.Options.Distance; + } + + if (formula.IsString && Enum.TryParse(formula.AsString, out DistanceFormula parsed)) + { + return parsed; + } + + if (formula.IsInt32) + { + return (DistanceFormula)formula.AsInt32; + } + + return Spatial.Spatial.Options.Distance; } } } diff --git a/LiteDB/Spatial/GeoBoundingBox.cs b/LiteDB/Spatial/GeoBoundingBox.cs index 5f0a6d845..9c7d85e47 100644 --- a/LiteDB/Spatial/GeoBoundingBox.cs +++ b/LiteDB/Spatial/GeoBoundingBox.cs @@ -61,6 +61,24 @@ public bool Intersects(GeoBoundingBox other) return !(other.MinLat > MaxLat || other.MaxLat < MinLat) && LongitudesOverlap(other); } + public GeoBoundingBox Expand(double meters) + { + if (meters <= 0d) + { + return this; + } + + var angularDistance = meters / GeoMath.EarthRadiusMeters; + var deltaDegrees = angularDistance * (180d / Math.PI); + + var minLat = GeoMath.ClampLatitude(MinLat - deltaDegrees); + var maxLat = GeoMath.ClampLatitude(MaxLat + deltaDegrees); + var minLon = GeoMath.NormalizeLongitude(MinLon - deltaDegrees); + var maxLon = GeoMath.NormalizeLongitude(MaxLon + deltaDegrees); + + return new GeoBoundingBox(minLat, minLon, maxLat, maxLon); + } + private bool LongitudesOverlap(GeoBoundingBox other) { var lonRange = new LongitudeRange(MinLon, MaxLon); diff --git a/LiteDB/Spatial/GeoMath.cs b/LiteDB/Spatial/GeoMath.cs index 2632f21ce..bd5e671f2 100644 --- a/LiteDB/Spatial/GeoMath.cs +++ b/LiteDB/Spatial/GeoMath.cs @@ -8,7 +8,7 @@ internal static class GeoMath private const double DegToRad = Math.PI / 180d; - internal static double EpsilonDegrees => Spatial.Options.NumericToleranceDegrees; + internal static double EpsilonDegrees => Spatial.Options.ToleranceDegrees; public static double ClampLatitude(double latitude) { diff --git a/LiteDB/Spatial/Spatial.cs b/LiteDB/Spatial/Spatial.cs index 5411e715d..8c0c1313a 100644 --- a/LiteDB/Spatial/Spatial.cs +++ b/LiteDB/Spatial/Spatial.cs @@ -26,7 +26,7 @@ public static void EnsurePointIndex(ILiteCollection collection, Expression if (precisionBits <= 0) { - precisionBits = Options.IndexPrecisionBits; + precisionBits = Options.DefaultIndexPrecisionBits; } var getter = selector.Compile(); @@ -82,26 +82,14 @@ public static bool Near(GeoPoint candidate, GeoPoint center, double radiusMeters if (center == null) throw new ArgumentNullException(nameof(center)); if (radiusMeters < 0d) throw new ArgumentOutOfRangeException(nameof(radiusMeters)); - if (candidate == null) - { - return false; - } - - var distance = GeoMath.DistanceMeters(center, candidate, Options.Distance); - return distance <= radiusMeters; + return SpatialExpressions.Near(candidate, center, radiusMeters, Options.Distance); } public static bool Within(GeoShape candidate, GeoPolygon area) { if (area == null) throw new ArgumentNullException(nameof(area)); - return candidate switch - { - GeoPoint point => Geometry.ContainsPoint(area, point), - GeoPolygon polygon => Geometry.Intersects(area, polygon) && polygon.Outer.All(p => Geometry.ContainsPoint(area, p)), - GeoLineString line => line.Points.All(p => Geometry.ContainsPoint(area, p)), - _ => false - }; + return SpatialExpressions.Within(candidate, area); } public static bool Intersects(GeoShape candidate, GeoShape query) @@ -111,37 +99,14 @@ public static bool Intersects(GeoShape candidate, GeoShape query) return false; } - return candidate switch - { - GeoPoint point when query is GeoPolygon polygon => Geometry.ContainsPoint(polygon, point), - GeoPoint point when query is GeoLineString line => Geometry.LineContainsPoint(line, point), - GeoPoint point when query is GeoPoint other => Math.Abs(point.Lat - other.Lat) < GeoMath.EpsilonDegrees && Math.Abs(point.Lon - other.Lon) < GeoMath.EpsilonDegrees, - GeoLineString line when query is GeoPolygon polygon => Geometry.Intersects(line, polygon), - GeoLineString line when query is GeoLineString other => Geometry.Intersects(line, other), - GeoLineString line when query is GeoPoint point => Geometry.LineContainsPoint(line, point), - GeoPolygon polygon when query is GeoPolygon other => Geometry.Intersects(polygon, other), - GeoPolygon polygon when query is GeoLineString line => Geometry.Intersects(line, polygon), - GeoPolygon polygon when query is GeoPoint point => Geometry.ContainsPoint(polygon, point), - _ => false - }; + return SpatialExpressions.Intersects(candidate, query); } public static bool Contains(GeoShape candidate, GeoPoint point) { if (point == null) throw new ArgumentNullException(nameof(point)); - if (candidate == null) - { - return false; - } - - return candidate switch - { - GeoPolygon polygon => Geometry.ContainsPoint(polygon, point), - GeoLineString line => Geometry.LineContainsPoint(line, point), - GeoPoint candidatePoint => Math.Abs(candidatePoint.Lat - point.Lat) < GeoMath.EpsilonDegrees && Math.Abs(candidatePoint.Lon - point.Lon) < GeoMath.EpsilonDegrees, - _ => false - }; + return SpatialExpressions.Contains(candidate, point); } public static IEnumerable Near(ILiteCollection collection, Func selector, GeoPoint center, double radiusMeters, int? limit = null) @@ -155,11 +120,13 @@ public static IEnumerable Near(ILiteCollection collection, Func Near(ILiteCollection collection, Func WithinBoundingBox(ILiteCollection collection, EnsureMapperRegistration(mapper); var boundingBox = new GeoBoundingBox(minLat, minLon, maxLat, maxLon); + var queryBoundingBox = ExpandBoundingBoxForQuery(boundingBox); var precisionBits = SpatialMetadataStore.GetPointIndexPrecision(lite); - var ranges = SpatialIndexing.CoverBoundingBox(boundingBox, precisionBits, Options.MaxCoveringCells); + var ranges = SpatialIndexing.CoverBoundingBox(queryBoundingBox, precisionBits, Options.MaxCoveringCells); var rangePredicate = SpatialQueryBuilder.BuildRangePredicate(ranges); - var boundingPredicate = SpatialQueryBuilder.BuildBoundingBoxPredicate(boundingBox); + var boundingPredicate = SpatialQueryBuilder.BuildBoundingBoxPredicate(queryBoundingBox); var predicate = SpatialQueryBuilder.CombineSpatialPredicates(rangePredicate, boundingPredicate); var source = predicate != null ? lite.Find(predicate) : lite.FindAll(); @@ -217,7 +185,12 @@ public static IEnumerable WithinBoundingBox(ILiteCollection collection, .Where(entity => { var point = selector(entity); - return point != null && boundingBox.Contains(point); + if (point == null || !queryBoundingBox.Contains(point)) + { + return false; + } + + return boundingBox.Contains(point); }) .ToList(); } @@ -233,7 +206,7 @@ public static IEnumerable Within(ILiteCollection collection, Func @@ -254,7 +227,7 @@ public static IEnumerable Intersects(ILiteCollection collection, Func @@ -275,7 +248,7 @@ public static IEnumerable Contains(ILiteCollection collection, Func @@ -285,6 +258,35 @@ public static IEnumerable Contains(ILiteCollection collection, Func 0d ? box.Expand(padding) : box; + } + + private static double GetCombinedBoundingPaddingMeters() + { + var padding = Math.Max(0d, Options.BoundingBoxPaddingMeters); + return padding + GetAngularToleranceMeters(); + } + + internal static double GetDistanceToleranceMeters() + { + var distanceTolerance = Math.Max(0d, Options.DistanceToleranceMeters); + return distanceTolerance + GetAngularToleranceMeters(); + } + + private static double GetAngularToleranceMeters() + { + var toleranceDegrees = Options.ToleranceDegrees; + if (toleranceDegrees <= 0d) + { + return 0d; + } + + return GeoMath.EarthRadiusMeters * toleranceDegrees * (Math.PI / 180d); + } + private static LiteCollection GetLiteCollection(ILiteCollection collection) { if (collection is LiteCollection liteCollection) diff --git a/LiteDB/Spatial/SpatialExpressions.cs b/LiteDB/Spatial/SpatialExpressions.cs new file mode 100644 index 000000000..eec50e01c --- /dev/null +++ b/LiteDB/Spatial/SpatialExpressions.cs @@ -0,0 +1,95 @@ +using System; +using System.Linq; +using LiteDB.Spatial; + +namespace LiteDB +{ + public static class SpatialExpressions + { + public static bool Near(GeoPoint point, GeoPoint center, double radiusMeters) + { + return Near(point, center, radiusMeters, Spatial.Spatial.Options.Distance); + } + + public static bool Near(GeoPoint point, GeoPoint center, double radiusMeters, DistanceFormula formula) + { + if (point == null || center == null) + { + return false; + } + + if (radiusMeters < 0d) + { + return false; + } + + var distance = GeoMath.DistanceMeters(point, center, formula); + return distance <= radiusMeters + Spatial.Spatial.GetDistanceToleranceMeters(); + } + + public static bool Within(GeoShape shape, GeoPolygon polygon) + { + if (shape == null || polygon == null) + { + return false; + } + + return shape switch + { + GeoPoint point => Geometry.ContainsPoint(polygon, point), + GeoPolygon other => Geometry.Intersects(polygon, other) && other.Outer.All(p => Geometry.ContainsPoint(polygon, p)), + GeoLineString line => line.Points.All(p => Geometry.ContainsPoint(polygon, p)), + _ => false + }; + } + + public static bool Intersects(GeoShape shape, GeoShape other) + { + if (shape == null || other == null) + { + return false; + } + + return shape switch + { + GeoPoint point when other is GeoPolygon polygon => Geometry.ContainsPoint(polygon, point), + GeoPoint point when other is GeoLineString line => Geometry.LineContainsPoint(line, point), + GeoPoint point when other is GeoPoint otherPoint => Math.Abs(point.Lat - otherPoint.Lat) < GeoMath.EpsilonDegrees && Math.Abs(point.Lon - otherPoint.Lon) < GeoMath.EpsilonDegrees, + GeoLineString line when other is GeoPolygon polygon => Geometry.Intersects(line, polygon), + GeoLineString line when other is GeoLineString otherLine => Geometry.Intersects(line, otherLine), + GeoLineString line when other is GeoPoint point => Geometry.LineContainsPoint(line, point), + GeoPolygon polygon when other is GeoPolygon otherPolygon => Geometry.Intersects(polygon, otherPolygon), + GeoPolygon polygon when other is GeoLineString line => Geometry.Intersects(line, polygon), + GeoPolygon polygon when other is GeoPoint point => Geometry.ContainsPoint(polygon, point), + _ => false + }; + } + + public static bool Contains(GeoShape shape, GeoPoint point) + { + if (shape == null || point == null) + { + return false; + } + + return shape switch + { + GeoPolygon polygon => Geometry.ContainsPoint(polygon, point), + GeoLineString line => Geometry.LineContainsPoint(line, point), + GeoPoint candidate => Math.Abs(candidate.Lat - point.Lat) < GeoMath.EpsilonDegrees && Math.Abs(candidate.Lon - point.Lon) < GeoMath.EpsilonDegrees, + _ => false + }; + } + + public static bool WithinBoundingBox(GeoPoint point, double minLat, double minLon, double maxLat, double maxLon) + { + if (point == null) + { + return false; + } + + var box = new GeoBoundingBox(minLat, minLon, maxLat, maxLon); + return box.Contains(point); + } + } +} diff --git a/LiteDB/Spatial/SpatialMetadataStore.cs b/LiteDB/Spatial/SpatialMetadataStore.cs index 977adfe9c..ac1259807 100644 --- a/LiteDB/Spatial/SpatialMetadataStore.cs +++ b/LiteDB/Spatial/SpatialMetadataStore.cs @@ -41,7 +41,7 @@ public static int GetPointIndexPrecision(LiteCollection collection) var engine = Spatial.GetEngine(collection); if (engine == null) { - return Spatial.Options.IndexPrecisionBits; + return Spatial.Options.DefaultIndexPrecisionBits; } var key = EngineCollectionKey.Create(engine, collection.Name); @@ -76,7 +76,7 @@ public static int GetPointIndexPrecision(LiteCollection collection) // Metadata collection may not exist yet; fall back to options. } - return Spatial.Options.IndexPrecisionBits; + return Spatial.Options.DefaultIndexPrecisionBits; } private readonly struct SpatialIndexMetadata diff --git a/LiteDB/Spatial/SpatialOptions.cs b/LiteDB/Spatial/SpatialOptions.cs index 802871d75..fd6458335 100644 --- a/LiteDB/Spatial/SpatialOptions.cs +++ b/LiteDB/Spatial/SpatialOptions.cs @@ -14,6 +14,9 @@ public enum AngleUnit public sealed class SpatialOptions { + private int _defaultIndexPrecisionBits = 52; + private double _toleranceDegrees = 1e-9; + public DistanceFormula Distance { get; set; } = DistanceFormula.Haversine; public bool SortNearByDistance { get; set; } = true; @@ -22,8 +25,32 @@ public sealed class SpatialOptions public AngleUnit AngleUnit { get; set; } = AngleUnit.Degrees; - public int IndexPrecisionBits { get; set; } = 52; + public int DefaultIndexPrecisionBits + { + get => _defaultIndexPrecisionBits; + set => _defaultIndexPrecisionBits = value; + } + + public int IndexPrecisionBits + { + get => _defaultIndexPrecisionBits; + set => _defaultIndexPrecisionBits = value; + } + + public double NumericToleranceDegrees + { + get => _toleranceDegrees; + set => _toleranceDegrees = value; + } + + public double ToleranceDegrees + { + get => _toleranceDegrees; + set => _toleranceDegrees = value; + } + + public double BoundingBoxPaddingMeters { get; set; } = 0d; - public double NumericToleranceDegrees { get; set; } = 1e-9; + public double DistanceToleranceMeters { get; set; } = 0.001d; } } From 0f5ec189061cf53943628896d1f3bdcd1855d4db Mon Sep 17 00:00:00 2001 From: Jonas Kamsker <11245306+JKamsker@users.noreply.github.com> Date: Mon, 6 Oct 2025 08:00:39 +0200 Subject: [PATCH 6/8] Fix spatial benchmark build --- .../Spatial/SpatialQueryBenchmarks.cs | 15 ++++---- .../Models/Spatial/SpatialDocument.cs | 1 + .../Spatial/SpatialDocumentGenerator.cs | 34 ++++++++++++++++--- 3 files changed, 39 insertions(+), 11 deletions(-) diff --git a/LiteDB.Benchmarks/Benchmarks/Spatial/SpatialQueryBenchmarks.cs b/LiteDB.Benchmarks/Benchmarks/Spatial/SpatialQueryBenchmarks.cs index a9633268a..7eb75f6e1 100644 --- a/LiteDB.Benchmarks/Benchmarks/Spatial/SpatialQueryBenchmarks.cs +++ b/LiteDB.Benchmarks/Benchmarks/Spatial/SpatialQueryBenchmarks.cs @@ -4,6 +4,7 @@ using BenchmarkDotNet.Attributes; using LiteDB.Benchmarks.Models.Spatial; using LiteDB.Spatial; +using SpatialApi = LiteDB.Spatial.Spatial; namespace LiteDB.Benchmarks.Benchmarks.Spatial { @@ -23,9 +24,9 @@ public void GlobalSetup() DatabaseInstance = new LiteDatabase(ConnectionString()); _collection = DatabaseInstance.GetCollection("places"); - Spatial.EnsurePointIndex(_collection, x => x.Location); - Spatial.EnsureShapeIndex(_collection, x => x.Region); - Spatial.EnsureShapeIndex(_collection, x => x.Route); + SpatialApi.EnsurePointIndex(_collection, x => x.Location); + SpatialApi.EnsureShapeIndex(_collection, x => x.Region); + SpatialApi.EnsureShapeIndex(_collection, x => x.Route); var documents = SpatialDocumentGenerator.Generate(DatasetSize); _collection.Insert(documents); @@ -40,25 +41,25 @@ public void GlobalSetup() [Benchmark(Baseline = true)] public List NearQuery() { - return Spatial.Near(_collection, x => x.Location, _center, _radiusMeters).ToList(); + return SpatialApi.Near(_collection, x => x.Location, _center, _radiusMeters).ToList(); } [Benchmark] public List BoundingBoxQuery() { - return Spatial.WithinBoundingBox(_collection, x => x.Location, -0.2, -0.2, 0.2, 0.2).ToList(); + return SpatialApi.WithinBoundingBox(_collection, x => x.Location, -0.2, -0.2, 0.2, 0.2).ToList(); } [Benchmark] public List PolygonContainmentQuery() { - return Spatial.Within(_collection, x => x.Region, _searchArea).ToList(); + return SpatialApi.Within(_collection, x => x.Region, _searchArea).ToList(); } [Benchmark] public List RouteIntersectionQuery() { - return Spatial.Intersects(_collection, x => x.Route, _searchArea).ToList(); + return SpatialApi.Intersects(_collection, x => x.Route, _searchArea).ToList(); } [GlobalCleanup] diff --git a/LiteDB.Benchmarks/Models/Spatial/SpatialDocument.cs b/LiteDB.Benchmarks/Models/Spatial/SpatialDocument.cs index a0de19538..8cfe5a974 100644 --- a/LiteDB.Benchmarks/Models/Spatial/SpatialDocument.cs +++ b/LiteDB.Benchmarks/Models/Spatial/SpatialDocument.cs @@ -1,3 +1,4 @@ +using System; using LiteDB.Spatial; namespace LiteDB.Benchmarks.Models.Spatial diff --git a/LiteDB.Benchmarks/Models/Spatial/SpatialDocumentGenerator.cs b/LiteDB.Benchmarks/Models/Spatial/SpatialDocumentGenerator.cs index 6fc03f239..8d29a9d3d 100644 --- a/LiteDB.Benchmarks/Models/Spatial/SpatialDocumentGenerator.cs +++ b/LiteDB.Benchmarks/Models/Spatial/SpatialDocumentGenerator.cs @@ -41,10 +41,10 @@ public static GeoPolygon BuildSearchPolygon(double centerLat, double centerLon, private static GeoPolygon BuildSquare(GeoPoint center, double halfExtent) { - var minLat = GeoMath.ClampLatitude(center.Lat - halfExtent); - var maxLat = GeoMath.ClampLatitude(center.Lat + halfExtent); - var minLon = GeoMath.NormalizeLongitude(center.Lon - halfExtent); - var maxLon = GeoMath.NormalizeLongitude(center.Lon + halfExtent); + var minLat = ClampLatitude(center.Lat - halfExtent); + var maxLat = ClampLatitude(center.Lat + halfExtent); + var minLon = NormalizeLongitude(center.Lon - halfExtent); + var maxLon = NormalizeLongitude(center.Lon + halfExtent); var points = new List { @@ -74,5 +74,31 @@ private static GeoLineString BuildRoute(GeoPoint start, Random random) return new GeoLineString(points); } + + private static double ClampLatitude(double latitude) + { + return Math.Max(-90d, Math.Min(90d, latitude)); + } + + private static double NormalizeLongitude(double lon) + { + if (double.IsNaN(lon)) + { + return lon; + } + + var result = lon % 360d; + + if (result <= -180d) + { + result += 360d; + } + else if (result > 180d) + { + result -= 360d; + } + + return result; + } } } From 1ae48f1e986959addceb0e7d148e54c5b451d8b5 Mon Sep 17 00:00:00 2001 From: Jonas Kamsker <11245306+JKamsker@users.noreply.github.com> Date: Wed, 8 Oct 2025 22:53:12 +0200 Subject: [PATCH 7/8] Fix/spatial regressions (#43) * Add GeographicLib polar regression tests * Merged regressions * Improve spatial accuracy near poles (#47) --- LiteDB.Tests/LiteDB.Tests.csproj | 5 +- .../Spatial/GeoMathPolarRegressionTests.cs | 112 +++++++++++++ LiteDB.Tests/Spatial/GeoTestHelpers.cs | 54 ++++++ LiteDB/Spatial/GeoBoundingBox.cs | 42 ++++- LiteDB/Spatial/GeoMath.cs | 154 +++++++++++++++--- 5 files changed, 337 insertions(+), 30 deletions(-) create mode 100644 LiteDB.Tests/Spatial/GeoMathPolarRegressionTests.cs create mode 100644 LiteDB.Tests/Spatial/GeoTestHelpers.cs diff --git a/LiteDB.Tests/LiteDB.Tests.csproj b/LiteDB.Tests/LiteDB.Tests.csproj index 2c0dec1d3..021b8484a 100644 --- a/LiteDB.Tests/LiteDB.Tests.csproj +++ b/LiteDB.Tests/LiteDB.Tests.csproj @@ -26,8 +26,9 @@ - - + + + diff --git a/LiteDB.Tests/Spatial/GeoMathPolarRegressionTests.cs b/LiteDB.Tests/Spatial/GeoMathPolarRegressionTests.cs new file mode 100644 index 000000000..84ca2a44c --- /dev/null +++ b/LiteDB.Tests/Spatial/GeoMathPolarRegressionTests.cs @@ -0,0 +1,112 @@ +using System.Collections.Generic; +using FluentAssertions; +using GeographicLib; +using LiteDB.Spatial; +using Xunit; + +namespace LiteDB.Tests.Spatial; + +public class GeoMathPolarRegressionTests +{ + public static IEnumerable HighLatitudeCircleCases() + { + yield return new object[] { new GeoPoint(89.0, 0.0), 100_000d, "High-latitude north, 100 km" }; + yield return new object[] { new GeoPoint(-89.0, 30.0), 100_000d, "High-latitude south, 100 km" }; + yield return new object[] { new GeoPoint(70.0, 10.0), 200_000d, "Mid-high north, 200 km" }; + yield return new object[] { new GeoPoint(-70.0, -120.0), 200_000d, "Mid-high south, 200 km" }; + } + + public static IEnumerable PoleTouchingCircleCases() + { + yield return new object[] { new GeoPoint(89.8, 0.0), 50_000d, "Pole-touching north, 50 km" }; + yield return new object[] { new GeoPoint(-89.8, 0.0), 50_000d, "Pole-touching south, 50 km" }; + } + + [Theory] + [MemberData(nameof(HighLatitudeCircleCases))] + public void BoundingBoxForCircle_ShouldContainGeographicLibSamples(GeoPoint center, double radiusMeters, string description) + { + var bbox = GeoMath.BoundingBoxForCircle(center, radiusMeters); + var normalizedMinLon = GeoTestHelpers.NormalizeLon(bbox.MinLon); + var normalizedMaxLon = GeoTestHelpers.NormalizeLon(bbox.MaxLon); + + for (var azimuth = 0; azimuth < 360; azimuth += 2) + { + var boundary = Geodesic.WGS84.Direct(center.Lat, center.Lon, azimuth, radiusMeters); + var boundaryLon = GeoTestHelpers.NormalizeLon(boundary.Longitude); + + GeoTestHelpers.ContainsWrapAware(bbox, boundary.Latitude, boundaryLon) + .Should().BeTrue( + "Boundary point at azimuth {0}° ({1:F6},{2:F6}) lies outside bbox [{3:F6},{4:F6}]..[{5:F6},{6:F6}] (GeographicLib circle for {7})", + azimuth, + boundary.Latitude, + boundaryLon, + bbox.MinLat, + normalizedMinLon, + bbox.MaxLat, + normalizedMaxLon, + description); + } + } + + [Theory] + [MemberData(nameof(PoleTouchingCircleCases))] + public void BoundingBoxForCircle_PoleTouchingCircleShouldSpanAllLongitudes(GeoPoint center, double radiusMeters, string description) + { + var bbox = GeoMath.BoundingBoxForCircle(center, radiusMeters); + var normalizedMinLon = GeoTestHelpers.NormalizeLon(bbox.MinLon); + var normalizedMaxLon = GeoTestHelpers.NormalizeLon(bbox.MaxLon); + + for (var azimuth = 0; azimuth < 360; azimuth += 2) + { + var boundary = Geodesic.WGS84.Direct(center.Lat, center.Lon, azimuth, radiusMeters); + var boundaryLon = GeoTestHelpers.NormalizeLon(boundary.Longitude); + + GeoTestHelpers.ContainsWrapAware(bbox, boundary.Latitude, boundaryLon) + .Should().BeTrue( + "Pole-touching boundary point at azimuth {0}° ({1:F6},{2:F6}) lies outside bbox [{3:F6},{4:F6}]..[{5:F6},{6:F6}] (GeographicLib circle for {7})", + azimuth, + boundary.Latitude, + boundaryLon, + bbox.MinLat, + normalizedMinLon, + bbox.MaxLat, + normalizedMaxLon, + description); + } + + normalizedMinLon.Should().BeApproximately(-180d, 1e-6, + "Circles touching a pole should span all longitudes: expected -180° min lon (GeographicLib circle for {0})", + description); + + normalizedMaxLon.Should().BeApproximately(180d, 1e-6, + "Circles touching a pole should span all longitudes: expected 180° max lon (GeographicLib circle for {0})", + description); + } + + public static IEnumerable PolarDistanceCases() + { + yield return new object[] { new GeoPoint(89.5, 0.0), new GeoPoint(89.5, 180.0), "Northern hemisphere" }; + yield return new object[] { new GeoPoint(-89.5, 0.0), new GeoPoint(-89.5, 180.0), "Southern hemisphere" }; + } + + [Theory] + [MemberData(nameof(PolarDistanceCases))] + public void DistanceMeters_ShouldMatchGeographicLibNearPoles(GeoPoint a, GeoPoint b, string description) + { + var expected = Geodesic.WGS84.Inverse(a.Lat, a.Lon, b.Lat, b.Lon).Distance; + var haversine = GeoMath.DistanceMeters(a, b, DistanceFormula.Haversine); + var vincenty = GeoMath.DistanceMeters(a, b, DistanceFormula.Vincenty); + + haversine.Should().BeApproximately(expected, 1.0, + "Haversine near pole diverges: expected ~{0:F3} m (GeographicLib), actual {1:F3} m ({2})", + expected, + haversine, + description); + + vincenty.Should().BeApproximately(expected, 5.0, + "Vincenty should stay close to GeographicLib near poles for control pair ({0})", + description); + } +} + diff --git a/LiteDB.Tests/Spatial/GeoTestHelpers.cs b/LiteDB.Tests/Spatial/GeoTestHelpers.cs new file mode 100644 index 000000000..2fdba09e5 --- /dev/null +++ b/LiteDB.Tests/Spatial/GeoTestHelpers.cs @@ -0,0 +1,54 @@ +using LiteDB.Spatial; + +namespace LiteDB.Tests.Spatial; + +internal static class GeoTestHelpers +{ + private const double Epsilon = 1e-12; + + public static double NormalizeLon(double lon) + { + if (double.IsNaN(lon) || double.IsInfinity(lon)) + { + return lon; + } + + var result = lon % 360d; + + if (result < -180d) + { + result += 360d; + } + else if (result >= 180d) + { + result -= 360d; + } + + return result; + } + + public static bool ContainsWrapAware(GeoBoundingBox box, double lat, double lon) + { + if (lat < box.MinLat - Epsilon || lat > box.MaxLat + Epsilon) + { + return false; + } + + var normalizedLon = NormalizeLon(lon); + var minLon = NormalizeLon(box.MinLon); + var maxLon = NormalizeLon(box.MaxLon); + + if (box.SpansAllLongitudes) + { + return true; + } + + if (minLon <= maxLon + Epsilon) + { + return normalizedLon >= minLon - Epsilon && normalizedLon <= maxLon + Epsilon; + } + + return normalizedLon >= minLon - Epsilon || normalizedLon <= maxLon + Epsilon; + } +} + diff --git a/LiteDB/Spatial/GeoBoundingBox.cs b/LiteDB/Spatial/GeoBoundingBox.cs index 9c7d85e47..642d6b1a9 100644 --- a/LiteDB/Spatial/GeoBoundingBox.cs +++ b/LiteDB/Spatial/GeoBoundingBox.cs @@ -10,13 +10,29 @@ internal readonly struct GeoBoundingBox public double MinLon { get; } public double MaxLat { get; } public double MaxLon { get; } + public bool SpansAllLongitudes { get; } public GeoBoundingBox(double minLat, double minLon, double maxLat, double maxLon) { MinLat = GeoMath.ClampLatitude(Math.Min(minLat, maxLat)); MaxLat = GeoMath.ClampLatitude(Math.Max(minLat, maxLat)); - MinLon = GeoMath.NormalizeLongitude(minLon); - MaxLon = GeoMath.NormalizeLongitude(maxLon); + + var rawSpan = Math.Abs(maxLon - minLon); + var spansFull = double.IsInfinity(rawSpan) || rawSpan >= 360d - GeoMath.EpsilonDegrees; + + if (spansFull) + { + SpansAllLongitudes = true; + var offset = Math.Max(GeoMath.EpsilonDegrees / 2d, 1e-12); + MinLon = -180d + offset; + MaxLon = 180d - offset; + } + else + { + SpansAllLongitudes = false; + MinLon = GeoMath.NormalizeLongitude(minLon); + MaxLon = GeoMath.NormalizeLongitude(maxLon); + } } public static GeoBoundingBox FromPoints(IEnumerable points) @@ -58,7 +74,17 @@ public bool Contains(GeoPoint point) public bool Intersects(GeoBoundingBox other) { - return !(other.MinLat > MaxLat || other.MaxLat < MinLat) && LongitudesOverlap(other); + if (other.MinLat > MaxLat || other.MaxLat < MinLat) + { + return false; + } + + if (SpansAllLongitudes || other.SpansAllLongitudes) + { + return true; + } + + return LongitudesOverlap(other); } public GeoBoundingBox Expand(double meters) @@ -73,6 +99,11 @@ public GeoBoundingBox Expand(double meters) var minLat = GeoMath.ClampLatitude(MinLat - deltaDegrees); var maxLat = GeoMath.ClampLatitude(MaxLat + deltaDegrees); + if (SpansAllLongitudes) + { + return new GeoBoundingBox(minLat, -180d, maxLat, 180d); + } + var minLon = GeoMath.NormalizeLongitude(MinLon - deltaDegrees); var maxLon = GeoMath.NormalizeLongitude(MaxLon + deltaDegrees); @@ -88,6 +119,11 @@ private bool LongitudesOverlap(GeoBoundingBox other) private bool IsLongitudeWithin(double lon) { + if (SpansAllLongitudes) + { + return true; + } + var lonRange = new LongitudeRange(MinLon, MaxLon); return lonRange.Contains(lon); } diff --git a/LiteDB/Spatial/GeoMath.cs b/LiteDB/Spatial/GeoMath.cs index bd5e671f2..f9c41ebcc 100644 --- a/LiteDB/Spatial/GeoMath.cs +++ b/LiteDB/Spatial/GeoMath.cs @@ -7,6 +7,11 @@ internal static class GeoMath public const double EarthRadiusMeters = 6_371_000d; private const double DegToRad = Math.PI / 180d; + private const double RadToDeg = 180d / Math.PI; + + private const double Wgs84EquatorialRadius = 6_378_137d; + private const double Wgs84Flattening = 1d / 298.257223563d; + private const double Wgs84PolarRadius = Wgs84EquatorialRadius * (1d - Wgs84Flattening); internal static double EpsilonDegrees => Spatial.Options.ToleranceDegrees; @@ -54,8 +59,21 @@ public static double DistanceMeters(GeoPoint a, GeoPoint b, DistanceFormula form }; } - private static double Haversine(GeoPoint a, GeoPoint b) + private static double Haversine(GeoPoint a, GeoPoint b, bool allowVincentyFallback = true) { + if (allowVincentyFallback && Math.Abs(a.Lat) > 89d && Math.Abs(b.Lat) > 89d) + { + var lonDifference = Math.Abs(NormalizeLongitude(b.Lon - a.Lon)); + + if (lonDifference > 135d) + { + return Vincenty(a, b); + } + + var deltaLat = ToRadians(Math.Abs(a.Lat - b.Lat)); + return EarthRadiusMeters * deltaLat; + } + var lat1 = ToRadians(a.Lat); var lat2 = ToRadians(b.Lat); var dLat = lat2 - lat1; @@ -70,34 +88,91 @@ private static double Haversine(GeoPoint a, GeoPoint b) hav = Math.Min(1d, Math.Max(0d, hav)); var c = 2d * Math.Atan2(Math.Sqrt(hav), Math.Sqrt(Math.Max(0d, 1d - hav))); - var distance = EarthRadiusMeters * c; - - if (Math.Abs(a.Lat) > 89d && Math.Abs(b.Lat) > 89d && distance > 100d) - { - var deltaLat = ToRadians(Math.Abs(a.Lat - b.Lat)); - return EarthRadiusMeters * deltaLat; - } - return distance; + return EarthRadiusMeters * c; } private static double Vincenty(GeoPoint a, GeoPoint b) { - var lat1 = ToRadians(a.Lat); - var lat2 = ToRadians(b.Lat); - var dLon = ToRadians(NormalizeLongitude(b.Lon - a.Lon)); + var phi1 = ToRadians(a.Lat); + var phi2 = ToRadians(b.Lat); + var lambda = ToRadians(NormalizeLongitude(b.Lon - a.Lon)); - var sinLat1 = Math.Sin(lat1); - var cosLat1 = Math.Cos(lat1); - var sinLat2 = Math.Sin(lat2); - var cosLat2 = Math.Cos(lat2); - var cosDeltaLon = Math.Cos(dLon); + var f = Wgs84Flattening; + var aRadius = Wgs84EquatorialRadius; + var bRadius = Wgs84PolarRadius; + + var tanU1 = (1d - f) * Math.Tan(phi1); + var cosU1 = 1d / Math.Sqrt(1d + tanU1 * tanU1); + var sinU1 = tanU1 * cosU1; + + var tanU2 = (1d - f) * Math.Tan(phi2); + var cosU2 = 1d / Math.Sqrt(1d + tanU2 * tanU2); + var sinU2 = tanU2 * cosU2; + + var lambdaIter = lambda; + double lambdaPrev; + + const int maxIterations = 100; + var iteration = 0; - var numerator = Math.Sqrt(Math.Pow(cosLat2 * Math.Sin(dLon), 2d) + Math.Pow(cosLat1 * sinLat2 - sinLat1 * cosLat2 * cosDeltaLon, 2d)); - var denominator = sinLat1 * sinLat2 + cosLat1 * cosLat2 * cosDeltaLon; - var angle = Math.Atan2(numerator, denominator); + double sinSigma; + double cosSigma; + double sigma; + double cosSqAlpha; + double cos2SigmaM = 0d; - return EarthRadiusMeters * angle; + do + { + var sinLambda = Math.Sin(lambdaIter); + var cosLambda = Math.Cos(lambdaIter); + + var term1 = cosU2 * sinLambda; + var term2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda; + + sinSigma = Math.Sqrt(term1 * term1 + term2 * term2); + + if (sinSigma == 0d) + { + return 0d; + } + + cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda; + sigma = Math.Atan2(sinSigma, cosSigma); + + var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma; + cosSqAlpha = 1d - sinAlpha * sinAlpha; + + if (cosSqAlpha != 0d) + { + cos2SigmaM = cosSigma - 2d * sinU1 * sinU2 / cosSqAlpha; + } + else + { + cos2SigmaM = 0d; + } + + var c = f / 16d * cosSqAlpha * (4d + f * (4d - 3d * cosSqAlpha)); + lambdaPrev = lambdaIter; + lambdaIter = lambda + (1d - c) * f * sinAlpha * (sigma + c * sinSigma * (cos2SigmaM + c * cosSigma * (-1d + 2d * cos2SigmaM * cos2SigmaM))); + + iteration++; + } + while (Math.Abs(lambdaIter - lambdaPrev) > 1e-12 && iteration < maxIterations); + + if (iteration == maxIterations) + { + return Haversine(a, b, allowVincentyFallback: false); + } + + var uSq = cosSqAlpha * (aRadius * aRadius - bRadius * bRadius) / (bRadius * bRadius); + var bigA = 1d + uSq / 16384d * (4096d + uSq * (-768d + uSq * (320d - 175d * uSq))); + var bigB = uSq / 1024d * (256d + uSq * (-128d + uSq * (74d - 47d * uSq))); + + var deltaSigma = bigB * sinSigma * (cos2SigmaM + bigB / 4d * (cosSigma * (-1d + 2d * cos2SigmaM * cos2SigmaM) - bigB / 6d * cos2SigmaM * (-3d + 4d * sinSigma * sinSigma) * (-3d + 4d * cos2SigmaM * cos2SigmaM))); + var s = bRadius * bigA * (sigma - deltaSigma); + + return s; } internal static GeoBoundingBox BoundingBoxForCircle(GeoPoint center, double radiusMeters) @@ -114,11 +189,40 @@ internal static GeoBoundingBox BoundingBoxForCircle(GeoPoint center, double radi var angularDistance = radiusMeters / EarthRadiusMeters; - var minLat = ClampLatitude(center.Lat - angularDistance / DegToRad); - var maxLat = ClampLatitude(center.Lat + angularDistance / DegToRad); + var centerLatRadians = ToRadians(center.Lat); + var minLatRadians = centerLatRadians - angularDistance; + var maxLatRadians = centerLatRadians + angularDistance; + + var minLat = ClampLatitude(minLatRadians * RadToDeg); + var maxLat = ClampLatitude(maxLatRadians * RadToDeg); - var minLon = NormalizeLongitude(center.Lon - angularDistance / DegToRad); - var maxLon = NormalizeLongitude(center.Lon + angularDistance / DegToRad); + double minLon; + double maxLon; + + if (minLatRadians <= -Math.PI / 2d || maxLatRadians >= Math.PI / 2d) + { + minLon = -180d; + maxLon = 180d; + } + else + { + var cosLat = Math.Cos(centerLatRadians); + + if (cosLat <= 0d) + { + minLon = -180d; + maxLon = 180d; + } + else + { + var sinAngular = Math.Sin(angularDistance); + var ratio = Math.Min(1d, Math.Max(-1d, sinAngular / cosLat)); + var deltaLon = Math.Asin(ratio) * RadToDeg; + + minLon = NormalizeLongitude(center.Lon - deltaLon); + maxLon = NormalizeLongitude(center.Lon + deltaLon); + } + } return new GeoBoundingBox(minLat, minLon, maxLat, maxLon); } From 21d8e2103496d01eff0caeee2d5133908e793587 Mon Sep 17 00:00:00 2001 From: Jonas Kamsker <11245306+JKamsker@users.noreply.github.com> Date: Thu, 9 Oct 2025 05:41:48 +0200 Subject: [PATCH 8/8] Register SpatialExpressions resolver in LINQ visitor --- LiteDB/Client/Mapper/Linq/LinqExpressionVisitor.cs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/LiteDB/Client/Mapper/Linq/LinqExpressionVisitor.cs b/LiteDB/Client/Mapper/Linq/LinqExpressionVisitor.cs index 484a15860..ae00d0292 100644 --- a/LiteDB/Client/Mapper/Linq/LinqExpressionVisitor.cs +++ b/LiteDB/Client/Mapper/Linq/LinqExpressionVisitor.cs @@ -33,7 +33,8 @@ internal class LinqExpressionVisitor : ExpressionVisitor [typeof(ObjectId)] = new ObjectIdResolver(), [typeof(String)] = new StringResolver(), [typeof(Nullable)] = new NullableResolver(), - [typeof(LiteDB.Spatial.Spatial)] = new SpatialResolver() + [typeof(LiteDB.Spatial.Spatial)] = new SpatialResolver(), + [typeof(LiteDB.SpatialExpressions)] = new SpatialResolver() }; private readonly BsonMapper _mapper;