From be25bca305cb30c14cc3b249833fbafeceb80479 Mon Sep 17 00:00:00 2001 From: Peter Stace Date: Sun, 1 Feb 2026 20:43:25 +1100 Subject: [PATCH] Optimize overlay operations with envelope checks and R-Tree indexing Add early exit for Intersection when envelopes are disjoint, returning an empty result of the appropriate dimension without calling the JTS port. Replace O(M*N) nested loops in gcAwareIntersection and gcAwareDifference with R-Tree indexed lookups for O(M log N) complexity. --- CHANGELOG.md | 7 ++++ geom/alg_overlay.go | 82 ++++++++++++++++++++++++++++++++++++++++----- geom/util.go | 9 ++++- 3 files changed, 89 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1155bbd6..fa44d2ec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,12 @@ # Changelog +## Unreleased + +- Optimize Intersection to return early when input envelopes are disjoint. + +- Optimize overlay operations (Intersection, Difference) for GeometryCollections + by using R-Tree indexing to reduce O(M×N) to O(M log N) complexity. + ## v0.57.0 2026-01-30 diff --git a/geom/alg_overlay.go b/geom/alg_overlay.go index c103a545..f2245750 100644 --- a/geom/alg_overlay.go +++ b/geom/alg_overlay.go @@ -1,6 +1,9 @@ package geom -import "github.com/peterstace/simplefeatures/internal/jtsport/jts" +import ( + "github.com/peterstace/simplefeatures/internal/jtsport/jts" + "github.com/peterstace/simplefeatures/rtree" +) // hasGC determines if either argument is a GeometryCollection. The JTS port // doesn't have full support for GeometryCollections, so we need to handle them @@ -16,6 +19,26 @@ func hasGC(a, b Geometry) bool { return a.IsGeometryCollection() || b.IsGeometryCollection() } +// envelopesDisjoint checks if two geometries have disjoint envelopes. +func envelopesDisjoint(a, b Geometry) bool { + return !a.Envelope().Intersects(b.Envelope()) +} + +// createEmptyResult creates an empty geometry with the given dimension. +// dim 0 = Point, dim 1 = LineString, dim 2 = Polygon, otherwise GeometryCollection. +func createEmptyResult(dim int) Geometry { + switch dim { + case 0: + return Point{}.AsGeometry() + case 1: + return LineString{}.AsGeometry() + case 2: + return Polygon{}.AsGeometry() + default: + return Geometry{} + } +} + // Union returns a geometry that represents the parts from either geometry A or // geometry B (or both). An error may be returned in pathological cases of // numerical degeneracy. @@ -39,6 +62,9 @@ func Intersection(a, b Geometry) (Geometry, error) { if a.IsEmpty() || b.IsEmpty() { return Geometry{}, nil } + if envelopesDisjoint(a, b) { + return createEmptyResult(minInt(a.Dimension(), b.Dimension())), nil + } if hasGC(a, b) { return gcAwareIntersection(a, b) } @@ -51,16 +77,33 @@ func gcAwareIntersection(a, b Geometry) (Geometry, error) { return Geometry{}, err } + // Build R-Tree from partsB. + items := make([]rtree.BulkItem, 0, len(partsB)) + for i, part := range partsB { + if box, ok := part.Envelope().AsBox(); ok { + items = append(items, rtree.BulkItem{Box: box, RecordID: i}) + } + } + tree := rtree.BulkLoad(items) + // The total result is the union of the intersections across the Cartesian // product of parts. var results []Geometry for _, partA := range partsA { - for _, partB := range partsB { - result, err := jtsOverlayOp(partA, partB, jts.OperationOverlayng_OverlayNG_INTERSECTION) + box, ok := partA.Envelope().AsBox() + if !ok { + continue + } + err := tree.RangeSearch(box, func(i int) error { + result, err := jtsOverlayOp(partA, partsB[i], jts.OperationOverlayng_OverlayNG_INTERSECTION) if err != nil { - return Geometry{}, err + return err } results = append(results, result) + return nil + }) + if err != nil { + return Geometry{}, err } } return UnaryUnion(NewGeometryCollection(results).AsGeometry()) @@ -118,20 +161,43 @@ func gcAwareDifference(a, b Geometry) (Geometry, error) { return Geometry{}, err } + // Build R-Tree from partsB. + items := make([]rtree.BulkItem, 0, len(partsB)) + for i, part := range partsB { + if box, ok := part.Envelope().AsBox(); ok { + items = append(items, rtree.BulkItem{Box: box, RecordID: i}) + } + } + tree := rtree.BulkLoad(items) + // The total result is the union of each part of A after each part of B has // been removed (sequentially). var results []Geometry for _, partA := range partsA { result := partA - for _, partB := range partsB { + box, ok := result.Envelope().AsBox() + if !ok { + continue + } + + err := tree.RangeSearch(box, func(i int) error { + // Recheck envelope intersection: result shrinks as differences are + // applied, so it may no longer overlap partsB[i]. + if envelopesDisjoint(result, partsB[i]) { + return nil + } var err error - result, err = jtsOverlayOp(result, partB, jts.OperationOverlayng_OverlayNG_DIFFERENCE) + result, err = jtsOverlayOp(result, partsB[i], jts.OperationOverlayng_OverlayNG_DIFFERENCE) if err != nil { - return Geometry{}, err + return err } if result.IsEmpty() { - break + return rtree.Stop } + return nil + }) + if err != nil { + return Geometry{}, err } results = append(results, result) } diff --git a/geom/util.go b/geom/util.go index 4ea5bc2d..0ed6bca5 100644 --- a/geom/util.go +++ b/geom/util.go @@ -7,7 +7,14 @@ import ( "sort" ) -// TODO: Remove this when we require Go 1.21 (the max builtin can be used instead). +// TODO: Remove these when we require Go 1.21 (the min/max builtins can be used instead). +func minInt(a, b int) int { + if a < b { + return a + } + return b +} + func maxInt(a, b int) int { if a > b { return a