Merge #62582

62582: geo/geomfn: add point in polygon optimization to st_covers r=otan a=andyyang890

This patch improves the performance of st_covers and
st_coveredby for the common use case of testing
whether a (multi)polygon covers a (multi)point.

Release note: None

Co-authored-by: Andy Yang <[email protected]>

pkg/geo/geomfn/binary_predicates.go

@@ -26,6 +26,19 @@ func Covers(a geo.Geometry, b geo.Geometry) (bool, error) {
 	if !a.CartesianBoundingBox().Covers(b.CartesianBoundingBox()) {
 		return false, nil
 	}
+
+	// Optimization for point in polygon calculations.
+	pointPolygonPair, pointKind, polygonKind := PointKindAndPolygonKind(a, b)
+	switch pointPolygonPair {
+	case PointAndPolygon:
+		// A point cannot cover a polygon.
+		return false, nil
+	case PolygonAndPoint:
+		// Computing whether a polygon covers a point is the same
+		// as computing whether a point is covered by the polygon.
+		return PointKindRelatesToPolygonKind(pointKind, polygonKind, PointPolygonCoveredBy)
+	}
+
 	return geos.Covers(a.EWKB(), b.EWKB())
 }
 
@@ -37,6 +50,17 @@ func CoveredBy(a geo.Geometry, b geo.Geometry) (bool, error) {
 	if !b.CartesianBoundingBox().Covers(a.CartesianBoundingBox()) {
 		return false, nil
 	}
+
+	// Optimization for point in polygon calculations.
+	pointPolygonPair, pointKind, polygonKind := PointKindAndPolygonKind(a, b)
+	switch pointPolygonPair {
+	case PolygonAndPoint:
+		// A polygon cannot be covered by a point.
+		return false, nil
+	case PointAndPolygon:
+		return PointKindRelatesToPolygonKind(pointKind, polygonKind, PointPolygonCoveredBy)
+	}
+
 	return geos.CoveredBy(a.EWKB(), b.EWKB())
 }
 
@@ -109,7 +133,8 @@ func Intersects(a geo.Geometry, b geo.Geometry) (bool, error) {
 
 	// Optimization for point in polygon calculations.
 	pointPolygonPair, pointKind, polygonKind := PointKindAndPolygonKind(a, b)
-	if pointPolygonPair {
+	switch pointPolygonPair {
+	case PointAndPolygon, PolygonAndPoint:
 		return PointKindRelatesToPolygonKind(pointKind, polygonKind, PointPolygonIntersects)
 	}
 
@@ -233,23 +258,41 @@ func Overlaps(a geo.Geometry, b geo.Geometry) (bool, error) {
 	return geos.Overlaps(a.EWKB(), b.EWKB())
 }
 
+// PointPolygonOrder represents the order of a point and a polygon
+// in an ordered pair of geometries.
+type PointPolygonOrder int
+
+const (
+	// NotPointAndPolygon signifies that a pair of geometries is
+	// not a point and a polygon.
+	NotPointAndPolygon PointPolygonOrder = iota
+	// PointAndPolygon signifies that the point appears first
+	// in an ordered pair of a point and a polygon.
+	PointAndPolygon
+	// PolygonAndPoint signifies that the polygon appears first
+	// in an ordered pair of a point and a polygon.
+	PolygonAndPoint
+)
+
 // PointKindAndPolygonKind returns whether a pair of geometries contains
 // a (multi)point and a (multi)polygon. It is used to determine if the
 // point in polygon optimization can be applied.
-func PointKindAndPolygonKind(a geo.Geometry, b geo.Geometry) (bool, geo.Geometry, geo.Geometry) {
-	var validPoint, validPolygon bool
-	var pointKind, polygonKind geo.Geometry
-	for _, geometry := range []geo.Geometry{a, b} {
-		switch geometry.ShapeType2D() {
-		case geopb.ShapeType_Point, geopb.ShapeType_MultiPoint:
-			validPoint = true
-			pointKind = geometry
+func PointKindAndPolygonKind(
+	a geo.Geometry, b geo.Geometry,
+) (PointPolygonOrder, geo.Geometry, geo.Geometry) {
+	switch a.ShapeType2D() {
+	case geopb.ShapeType_Point, geopb.ShapeType_MultiPoint:
+		switch b.ShapeType2D() {
 		case geopb.ShapeType_Polygon, geopb.ShapeType_MultiPolygon:
-			validPolygon = true
-			polygonKind = geometry
+			return PointAndPolygon, a, b
+		}
+	case geopb.ShapeType_Polygon, geopb.ShapeType_MultiPolygon:
+		switch b.ShapeType2D() {
+		case geopb.ShapeType_Point, geopb.ShapeType_MultiPoint:
+			return PolygonAndPoint, b, a
 		}
 	}
-	return validPoint && validPolygon, pointKind, polygonKind
+	return NotPointAndPolygon, a, b
 }
 
 // Touches returns whether geometry A touches geometry B.

pkg/geo/geomfn/binary_predicates_test.go

@@ -42,6 +42,7 @@ var (
 	leftRectHoleCornerPoint = geo.MustParseGeometry("POINT(-0.75 0.75)")
 	leftRectHoleEdgePoint   = geo.MustParseGeometry("POINT(-0.75 0.5)")
 	leftRectMultiPoint      = geo.MustParseGeometry("MULTIPOINT(-0.5 0.5, -0.9 0.1)")
+	leftRectEdgeMultiPoint  = geo.MustParseGeometry("MULTIPOINT(-1.0 0.2, -0.9 0.1)")
 )
 
 func TestCovers(t *testing.T) {
@@ -53,6 +54,9 @@ func TestCovers(t *testing.T) {
 		{rightRect, rightRectPoint, true},
 		{rightRectPoint, rightRect, false},
 		{leftRect, rightRect, false},
+		{leftRect, leftRectEdgePoint, true},
+		{leftRectWithHole, leftRectPoint, false},
+		{leftRect, emptyPoint, false},
 	}
 
 	for i, tc := range testCases {
@@ -78,6 +82,8 @@ func TestCoveredBy(t *testing.T) {
 		{rightRect, rightRectPoint, false},
 		{rightRectPoint, rightRect, true},
 		{leftRect, rightRect, false},
+		{leftRectEdgeMultiPoint, leftRectWithHole, true},
+		{leftRectPoint, emptyRect, false},
 	}
 
 	for i, tc := range testCases {
@@ -110,6 +116,8 @@ func TestContains(t *testing.T) {
 		{leftRectWithHole, leftRectMultiPoint, false},
 		{leftRect, leftRectMultiPoint, true},
 		{bothLeftRectsHoleFirst, leftRectPoint, true},
+		{leftRect, leftRectEdgePoint, false},
+		{leftRect, leftRectEdgeMultiPoint, true},
 	}
 
 	for i, tc := range testCases {

pkg/geo/geomfn/distance.go

@@ -844,6 +844,9 @@ const (
 	// PointPolygonIntersects is the relationship where a (multi)point
 	// intersects a (multi)polygon.
 	PointPolygonIntersects PointPolygonRelationType = iota + 1
+	// PointPolygonCoveredBy is the relationship where a (multi)point
+	// is covered by a (multi)polygon.
+	PointPolygonCoveredBy
 )
 
 // PointKindRelatesToPolygonKind returns whether a (multi)point and
@@ -864,6 +867,9 @@ func PointKindRelatesToPolygonKind(
 
 	// Check whether each point intersects with at least one polygon.
 	// For Intersects, at least one point must intersect with at least one polygon.
+	// For CoveredBy, every point must intersect with at least one polygon.
+	intersectsOnce := false
+pointOuterLoop:
 	for {
 		point, hasPoint, err := pointKindIterator.Next()
 		if err != nil {
@@ -887,16 +893,29 @@ func PointKindRelatesToPolygonKind(
 				return false, err
 			}
 			if intersects {
+				intersectsOnce = true
 				switch relationType {
 				case PointPolygonIntersects:
+					// A single intersection is sufficient.
 					return true, nil
+				case PointPolygonCoveredBy:
+					// If the current point intersects, check the next point.
+					continue pointOuterLoop
 				default:
 					return false, errors.Newf("unknown PointPolygonRelationType")
 				}
 			}
 		}
+		// Case where a point in the (multi)point does not intersect
+		// a polygon in the (multi)polygon.
+		switch relationType {
+		case PointPolygonCoveredBy:
+			// Each point in a (multi)point must intersect a polygon in the
+			// (multi)point to be covered by it.
+			return false, nil
+		}
 	}
-	return false, nil
+	return intersectsOnce, nil
 }
 
 // pointIntersectsPolygon returns whether a point intersects with a polygon.
@@ -915,7 +934,7 @@ func pointIntersectsPolygon(point geom.T, polygon geom.T) (bool, error) {
 		return false, errors.Newf("geomToGeodist failed to convert a *geom.Point to a *geodist.Point")
 	}
 
-	// Convert polygon from a geo.Geometry to a geodist.Polygon.
+	// Convert polygon from a geom.T to a geodist.Polygon.
 	_, ok = polygon.(*geom.Polygon)
 	if !ok {
 		return false, errors.Newf("second geometry passed to PointIntersectsPolygon must be a polygon")