Skip to content

EdgeQuery with MinDistanceToShapeIndexTarget fails to find shapes within a distance limit when using optimized method #293

Description

@joshua-zingale-rincon

Issue

A query constructed with NewClosestEdgeQuery against a MinDistanceToShapeIndexTarget target fails to find non-overlapping shapes within the limit distance when the optimized (non-brute-force) method is used.

I am not able to commit time to fixing this at the moment, though I may give a quick stab at it.

Issue Context & Implications

I was using EdgeQuery to find all shapes within a distance of a given shape to determine neighbors in a clustering algorithm.
This issue caused significant and pervasive silent failures to merge nearby shapes into the same cluster.

Failing Test Code

package s2edgequeryfailure

import (
	"math"
	"testing"

	"github.com/golang/geo/s1"
	"github.com/golang/geo/s2"
)

// Bug: with MinDistanceToShapeIndexTarget, EdgeQuery's optimized
// path (the default) fails to find a neighbor that the brute-force path finds,
// for the same index, target, and distance limit.
//
// shapeA and shapeB are two small quads < 20 m apart. A query with a 20 m
// DistanceLimit should report shapeB as within range of shapeA. It does under
// UseBruteForce(true), but returns nothing under the default optimized path.
//
// The optimized path is only taken once the index holds more edges than the
// target's brute-force threshold, so multiple far-away filler quads are added to
// push the index over it. They sit >100 km away and never match.
func TestShapeIndexTargetOptimizedPathMissesNeighbor(t *testing.T) {
	// Two ~2 m quads whose nearest edges are ~9 m apart, north-south.
	// Anchored at the real coordinates that reproduced the bug.
	const (
		quadHalf = 1.0 // 2 m squares
		gap      = 9.0 // ~9 m edge-to-edge, well within the 20 m limit
	)
	baseLat, baseLng := 14.56145, -90.53163
	centerOffset := quadHalf + gap/2 // center-to-center = 2*half + gap

	// The target shape
	shapeA := quadAround(baseLat-metersToLatDegrees(centerOffset), baseLng, quadHalf, quadHalf)
	shapeB := quadAround(baseLat+metersToLatDegrees(centerOffset), baseLng, quadHalf, quadHalf)

	// Far-away filler (>100 km north), just to raise the index edge count so
	// the query takes the optimized path instead of brute force.
	index := s2.NewShapeIndex()
	for i := 1; i < 7; i++ {
		index.Add(quadAround(baseLat+float64(i), baseLng, 1, 1))
	}
	idB := index.Add(shapeB)
	index.Build()

	// Target: a one-shape index holding shapeA
	targetIndex := s2.NewShapeIndex()
	targetIndex.Add(shapeA)
	targetIndex.Build()
	target := s2.NewMinDistanceToShapeIndexTarget(targetIndex)

	limit := s1.ChordAngleFromAngle(s1.Angle(20.0 / 6371000.0)) // 20 m

	find := func(useBruteForce bool) bool {
		q := s2.NewClosestEdgeQuery(index, s2.NewClosestEdgeQueryOptions().
			DistanceLimit(limit).
			IncludeInteriors(true).
			UseBruteForce(useBruteForce))
		for _, r := range q.FindEdges(target) {
			if r.ShapeID() == idB {
				return true
			}
		}
		return false
	}

	brute := find(true)
	optimized := find(false)

	if !brute {
		t.Fatal("brute-force path did not find shapeB; fixture broken")
	}
	if !optimized {
		t.Errorf("optimized path did NOT find shapeB within 20 m, but brute force did "+
			"(shapeB id=%d) — EdgeQuery optimized/brute-force mismatch with "+
			"MinDistanceToShapeIndexTarget", idB)
	}
}

const earthRadiusMeters = 6371000.0

// metersToLatDegrees converts a north-south distance in meters to degrees of
// latitude
func metersToLatDegrees(m float64) float64 {
	return m / (earthRadiusMeters * math.Pi / 180.0)
}

// metersToLngDegrees converts an east-west distance in meters to degrees of
// longitude at the given latitude.
func metersToLngDegrees(m, atLatDeg float64) float64 {
	return m / (earthRadiusMeters * math.Pi / 180.0 * math.Cos(atLatDeg*math.Pi/180.0))
}

// quadAround builds a CCW axis-aligned quad centered at (latDeg, lngDeg)
// with the given half-width and half-height in meters.
func quadAround(latDeg, lngDeg, halfW, halfH float64) *s2.Polygon {
	dLat := metersToLatDegrees(halfH)
	dLng := metersToLngDegrees(halfW, latDeg)
	latLo, latHi := latDeg-dLat, latDeg+dLat
	lngLo, lngHi := lngDeg-dLng, lngDeg+dLng
	return s2.PolygonFromLoops([]*s2.Loop{s2.LoopFromPoints([]s2.Point{
		s2.PointFromLatLng(s2.LatLngFromDegrees(latLo, lngLo)),
		s2.PointFromLatLng(s2.LatLngFromDegrees(latLo, lngHi)),
		s2.PointFromLatLng(s2.LatLngFromDegrees(latHi, lngHi)),
		s2.PointFromLatLng(s2.LatLngFromDegrees(latHi, lngLo)),
	})})
}

Go Mod

go 1.26.2

require github.com/golang/geo v0.0.0-20260713102120-857a528af641

require github.com/google/go-units v0.0.0-20250612230646-eddd77f68220 // indirect

AI Use

I found and reproduced this bug by hand but used Claude Opus 4.8 to generate the final test case above.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions